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ABSTRACT 


The dynamic stall of rapidly pitching and oscillating 
airfoils is investigated by the numerical solution of the 
full compressible unsteady Lsv-dimensionali Navier-stoKxes 
equations using an alternating-direction-implicit scheme. 
The flow is assumed to be fully turbulent, and the turbulent 
stresses are modelled by the Baldwin-Lomax eddy viscosity 
mode]. Three airfoils (NACA 0012, NACA 0012-33, and NACA 
0012-63) are analyzed for the purpose of examining the 
influence of leading-edge geometry on unsteady flow separa- 
tion. It is found that a larger leading edge radius, thick- 
er contouring of the forward part of the airfoil, or in- 
creasing reduced frequency results in delaying flow separa- 
tion and formation of the dynamic stall vortex to a higher 
angle of attack, yielding higher peak C.. Within the scope 
of this study, the pressure gradient encountered by the flow 
at initial separation is found to be independent of reduced 
frequency and freestream speed. The critical pressure 
gradient is dependent on leading edge radius and increases 


for decreasing leading edge radius. 
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I. INTRODUCTION 


A. BACKGROUND 

The aeronautical community has been aware of the phenome- 
non of dynamic stall for several decades. Dynamic stall is 
characterized b, separated flow and shedding of the leading 
edge vortex from the upper surface of an airfoil which is 
rapidly pitched to angles of attack beyond the normal stalling 
angle of attack. The phenomenon results in significant 
temporary increases and decreases in lift, drag, and moment 
coefficients. 

The scope of this work is the investigation of the 
phenomenon of dynamic stall for rapidly pitching and oscillat- 
ing airfoils. The phenomenon of dynamic stall was observed 
for the first time for flows over the retreating blades of a 
helicopter. The dynamic stall resulting from the oscillatory 
motion of the rotor blade is associated with an increase in 
lift and the development of a severe nose down pitching 
momer.t. The effects of dynamic stall are usually undesirable 
for the helicopter, and where possible special care is taken 
to reduce its effects by special design of the rotor. On the 
other hand, interest has recently developed to exploit the 
increased lift obtained from the rapid pitch-up motion of an 


airfoil in order to enhance the maneuverability and extend the 





flight envelope of the modern fighter aircraft to the high 
angle of attack regime, or to alleviate retreating blade stall 
for helicopters utilizing Higher Harmonic Control (HHC). 

It is well known that for flow over an airfoil at fixed 
angle of attack the streamlined airflow is disrupted once a 
critical angle of attack is exceeded. At stall, the flow over 
the upper surface of the airfoil separates and the lift drops. 
It was observed, however, that rapid pitch-up motion of the 
airfoil delays static stall so that high lift .an be main- 
tained for angles of attack beyond the static stall angle. 

As the pitch-up angle exceeds the static stalling angle of 
attack, a thin layer of reversed flow develops in the boundary 
layer. This reversed flow occurs in two types of stall: 
trailing edge stall and leading edge stall [Ref. 1]. In 
trailing edge stall, the reverse flow region begins near the 
trailing edge and traverses forward; in leading edge stall, 
the reverse flow region occurs near the suction peak just aft 
of the leading edge. In both cases, a vortex begins to form 
near the leading edge region, expands, and moves downstream. 
The angle of attack at which the vortex is formed depends on 
airfoil shape, pitch rate, mean angle and amplitude, Mach 
number, and Reynolds number. 

The vortex formed at the leading edge is called the 
dynamic stall vortex, and moves with a speed of approximately 
0.4 freestream speed relative to the airfoil as the pitch-up 


progresses. Lift, drag, and pitching moment increase signifi- 





cantly until the vortex approaches the trailing edge, then 
drop sharply, but not simultaneously (Figure 1). The unsteady 
surface pressure increases, and the suction peak appears at 
different locations along the chord as the dynamic stall 
vortex moves over the airfoil. Secondary and tertiary 
vortices may also be present and produce additional suction 
peaks and fluctuations in airloads. 


To date, most theo- 
(a) STATIC STALL ANGLE EXCEEDED 


{b) FIRST APPEARANCE OF FLOW 
retical investigations REVERSAL ON SURFACE 
4 
Sa 
ic) LARGE EQOIES APPEAR IN 
BOUNDARY LAYER 


focused on the effects NG 


. : id) FLOW REVERSAL SPREADS OVER 
of variations of reduced \ MUCH OF AIRFOIL CHORD 


frequency (k = &c/2U,), te) Se 
\ 


on dynamic stall have 


{eo} VORTEX FORMS NEAR 
LEADING EOGE 


SS 


angle of attack, free 


z 
v 
w 
oO 
« 
°o 
u 
2 
< 
= 
4 
° 
z 


stream speed, and 


Reynolds number on a 
{g)} MOMENT STALL OCCURS 


(h) LIFT STALL BEGINS 
{i) MAXIMUM NEGATIVE MOMENT 


particular airfoil shape 
(usually the NACA 0012 


airfoil). McCroskey 


PITCHING MOMENT, Cy 


{[Ref. 1] documented the 


(kh) BOUNDARY LAYER REATTACHES, 
FRONT TO REAR 


effects of reduced fre- EEE aE = 


INCIDENCE, a, deg fl) RETUAN TO UNSTALLED VALUES 





quency, amplitude and 
Figure 1. Events of Dynamic Stall 

Mach number on different on NACA 0012 Airfoil 

airfoils. These experi- 


mental studies demonstrated the large hysteresis in the 





occurrence of stall of an oscillating airfoil compared to a 
fixed angle of attack airfoil. McCroskey and Pucci (Ref. 2] 
identified varying regimes of viscous-inviscid interacticn 
during varying degrees of unsteady flow separation. The 
conclusion of the experimental studies of Refs. 1 and 2 was 
that the reduced frequency has a dominant effect on the 
development and progression of the dynamic stall. Experi- 
mental work by Chandrasekhara and Carr [{Ref. 3] using the NACA 
0012 airfoil showed that a dynamic stall vortex always forms 
near the leading edge of an oscillating airfoil. Their study 
also documents the movement of the dynamic stall vortex. 
Chandrasekhara and Brydges [Ref. 4] documented the effects of 
increasing amplitude on an oscillating airfoil in both 
compressible and incompressible flow. They showed that larger 
amplitudes resulted in vortex retention at higher angles of 
attack for a given Mach number and reduced freguency. 
Progress in computational fluid dynamics has made possible 
the study of dynamic stall by numerical solution of the 
unsteady Navier-Stokes equations. Mehta [Ref. 5] demonstrated 
that Navier-Stokes simulation of the unsteady incompressible 
flow around the airfoil in oscillatory motion can reproduce 
the experimentally observed results. Wu et al. [Ref. 6] 
presented solution procedures based on an integral formulation 
of the incompressible Navier-Stokes equations for the 
computation of unsteady flow over airfoils. Beddoes [Ref. 7] 


and Jang et al. [Ref. 8] presented viscous-inviscid computa- 








tion methods for unsteady flows. Sankar and Tang [Ref. 9], 
Visbal [Ref. 10], and Ekaterinaris [Ref. 11] used ADI (Alter- 
nating Direction Implicit) numerical schemes for the solution 
of the compressible Navier-Stokes equations, and investigated 
the effects of compres ibility on dynamic stall. The numeri- 
cal solutions for both incompressible and compressible flows 
showed good agreement with experimental results, and predicted 


the events of dynamic stalls. 


B. PURPOSE 

As indicated, the aforementioned studies focused on the 
variation of flowfield parameters on the dynamic stall of a 
particular airfoil. The purpose of this study is the system- 
atic investigation of the effect of leading edge geometry on 
the development of dynamic stall. It is expected that the 
effect of the leading edge geometry will have a significant 
effect on both the development of the dynamic stall vortex and 
its subsequent shedding. This numerical investigation is 
intended to provide a cost-effective means of quantifying 
which airfoil parameter variations should be analyzed in more 
expensive wind-tunnel tests. The numerical investigation of 
airfoil parameter variations offers the benefit of optimizing 
the utilization of more costly test facilities. 

The investigation is conducted using a numerical solution 


of the full two-dimensional, Reynolds-averaged Navier-Stokes 








Equations, and the Baldwin-Lomax eddy viscosity model of Ref. 


16 is used to obtain the turbulent stresses. 


C. AIRFOIL SELECTION 

Modifications to the NACA 0012 airfoil were made forward 
of the point of maximum thickness (12% thick at 30%) chord by 
the method suggested in Ref. 12. The NACA 0012 profile was 
retained aft of 30% chord for the two new airfoils. In this 
manner, changes in the flow could be attributed directly to 
the changes in the forward part of the airfoil. The two air- 
foils consisted of the NACA 0012-63 and NACA 0012-33 airfoil 
section forward of 30% chord. The NACA 0012-63 has the same 
leading edge radius as the NACA 0012 (1.58% chord) but 
different curvature to the point of maximum thickness. The 
NACA 0012-33 has a smaller leading edge radius (0.39% chord) 
with curvature necessary to achieve the same 12% thickness 
ratio at 30% chord. Airfoil coordinates are provided in Table 


1, and the resulting modified profiles are shown in Figure 2. 











TABLE 1. AIRFOIL COORDINATES 
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-O1 -01701 
-02 - 02360 
-03 -02842 
-04 .03232 
-05 -03555 
- 06 -03840 
.O7 -04090 
.08 -04313 
-09 -04505 
-10 -04681 
-li1 .04842 
~12 -04990 
.13 -05123 
.14 05242 
-15 .05345 
-16 -05441 
.18 -05610 
20 -05742 
-22 -05840 
.24 .05903 
.25 -05943 
26 -05964 
28 -05992 
.30 -06000 
-40 -05803 
-50 .05294 
- 60 .04563 
- 70 -03664 
80 -02623 
-90 01448 
-95 -00807 
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NO012-63 


.01682 
.02320 
.02783 
.03160 
.03474 
.03750 
.03991 
.04210 
.04403 
.04580 
.04742 
.04890 
.05023 
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.00807 
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Airfoil Profile Comparison 


Figure 2. 





II. GOVERNING EQUATIONS 


The flow of a compressible, viscous fluid satisfies 
conservation of mass, momentum, and energy. The conservation 
of mass is expressed by the Continuity Equation, the conserva- 
tion of momentum by the Navier-Stokes Equations, and the 
conservation of energy by the Energy Equation. The conserva- 
tion equations are derived in the following section. The 


derivation is done in an Eulerian frame of reference. 


A. CONTINUITY EQUATION 

Consider a control volume in a flow field where flow 
properties vary with both time and space. Conservation of 
mass requires that the rate of change of mass inside the 
control volume equals net mass flux out of the control volume. 


This is expressed as: 


ZSff.e d(vol) + [ Jovas = 0 (1) 


using Gauss's theorem for a surface integral, Eq. 1 becomes 


(ff. acvoa) + fff Ve) a(vo1) = 0 (2) 


or, 


fff, '2. + V(pWld(vol) = 0 (3) 


therefore, for an arbitrary control volume, 








OP esq, 7 
ae * V(pY) 0 (4) 


for any continuous flow. For a two-dimensional flow, and 


Cartesian coordinates, Eq. 4 reduces to 


Go , A(pu) r O(pv) _ 0 (5) 
ot Ox oy 


B. THE MOMENTUM (NAVIER-STOKES) EQUATIONS 

Consider a fluid element in a rectangular Cartesian 
coordinate system. The stresses and pressures are shown in 
Figure 3, and body forces are neglected. Summation of forces 


in the x-direction yields: 


Du 
F,= 9 dx dy dz—- 

















 iPpetpacull OF ea (6) 
= [p-(p et clea [(t 4 FF ee Tx) Gydz 
+ [(tyy+ FEA) ~ Fy dxdz + [t.,+ “2% dz) -1,,] dxdy 
or, 
po oye ae OF sxx + OF yx + OY ax (7) 
Dt Ox Ox oy Oz 


Similarly, summation of forces in the y- and z- directions 








yields: 
Dy. op, Mw, Ay , Ky (8) 
PDE oy "ox 7 oy * “Gz 
Dw. Op , Me , Nye , Nee (9) 
P Dt dz Ox oy * “oz 
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components 


1,, dx dy | 





p+ 2 as) ) ay dz 


ete + was ax) dy dz 








Figure 3. Forces on a Fluid Elemeni 


Using the Continuity Equation, the equations expressing the 


conservation of momentum for a two-dimensional flow are: 








a a 
Seu . 2 (pu? + p) + £ (pur) rat, Sr (10) 





otpy) 4 ZS (puv) + (pv? + p) cae ot (11) 


and f¢ to the 


The relation of the viscous stresses ft yy’ i 


rr ae | 
independent variables is developed in the following discus- 
Sion. 


For a two-dimensional flow in Cartesian coordinates with 


an infinitesimal fluid element undergoing distortion due to 
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Fluid clement at 
tune. ¢ + Ot 


\ 
Fluid element at time ¢ 





Figure 4. Strain ona Fluid Element 


stresses as shown in Figure 4, the angular displacements A®@, 


and A@, are: 


= 40Y = -9u 12 
A8, = SY At AO, = -sAt (12) 

The strain increment is given by: 
AK = A6, - AB, , (13) 


in the limit, the rate of strain is given by: 
ge vias RR nay Bod ED tel Ng 2 (14) 


Using Newton's Law of Fluid Friction (definition of 


viscosity), 


Ty 2 ba (15) 

















then the rate of strain caused by the tangential shear stress 


is given by: 


= = ponte Riatind 16 
Ty = BE gy Wa + ay! (16) 


For large velocity gradients the normal stresses t, and Ty 
can be sianificant and result in a viscous-induced normal 
force on the fluid element. For example, as documented by 
Schlichting [Ref. 13], in order for fluid isotropy to be 
maintained at every point, the principal axes of stress and 
rate-of-strain must coincide to avoid introducing a preferred 
rotation direction. With this concept, a normal stress must 
depend both on its respective component rate of strain as well 
as the shearing strain rates, with different weighting 
factors. Choosing the factor 2p for the component direction 


factor, which causes Newton's Law of Friction to be satisfied 


for simple shear, obtain 


i Ou. Ov Ou 17 
Tux Riss ho + A Fe ( ) 
Ou. Ov Ov 
= a Pnclialt 18 
© ye A ax ae + 2p ay (18) 


where 4 is the shearing stress proportionality factor. Using 


Stokes's hypothesis, 4 = -2/3 yj, obtain 
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a 
u 
{ 


ar 3 ox oy 
Ty = 4 (28¥ _ Su) (19) 
Ov, du 
T xy p ( oi oF 


C. THE ENERGY EQUATION 

Conservation of energy is the manifestation of the First 
Law of Thermodynamics (dE = dW + dQ) to a moving fluid element 
in rectangular Cartesian coordinates. The First Law of 
Thermodynamics is applied to a control volume (Figure 5), 
where the energy fluxes are shown. The rate of change of 
energy inside the fluid element is equal to the net flux of 
heat into the element plus the rate of work done on the 
element by pressure and viscous stress. 

The rate of change of energy of the fluid element having 
an instantaneous internal energy per unit mass e and speed V 


is 


poles 1) axdydz (20) 


where Vi =u'+v'+w. The heat flux into the fluid element 
is the sum of external volumetric heating and heat transfer 
across the surface due to temperature gradients. Assuming no 
external heat addition, the volumetric heating is zero. 

The net heat transferred out of the fluid element due to 


thermal conduction in the x-direction can be expressed as: 
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[up + Buu) ds} dy dz 


O(ur,.) 
(urs + a ae ds] dy dz 
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oq 
(. + a as) dy dz 


a(ur,, ) P 


Bz 7] dx dy 





Figure 5. Energy Fluxes on Fluid Element 


SS ts 0d, _ _ 9d, i 
Lgy= Ag4 aye ax) ] dydz = 5 dxdydz (21) 








Accounting for the y and z directions, the total heat trans- 


ferred out of the fluid element is: 


oq, og, + od, 


“Oy = AWE 22 
ax | oy ag ) exdydz (V-q) dxdy dz (22) 


mA 








where the heat flux is expressed in terms of the temperature 


gradient according to Fourier's law of heat conduction as 


_, or (23) 
dy = ko 


and k is the thermal conductivity, considered constant and 
independent of temperature. 
The rate of work done on the fluid element in the 


x-direction by pressure and shear stresses is 
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_O(up) , 0 (ut ,,) O(ut ,,) 0 (ut ,,) a 
ogee se axce ey  gg eee 


Accounting for the y and z direction pressure and shear 
stresses, the net rate of work done on the fluid element can 


be expressed as: 


z, DUT) | Aluty,) | Aut zy) 


IVipy Ox oy oz 
, OMT ay) , OLY yy) , GAM zy) (25) 
ox Oy oz 
re) ro] rs) 


The complete energy equation for a viscous flow with no 


external heating can then be expressed as: 


p— (e+ Sv?) = [-V-q-VipvV 
. A(ut ,,) P d(ut,,) : d(ut ,,) 


Ox oy Oz P 
O(vt,,) . A(vty,) | A(vt zy) oie 
+ — + St 
Ox oy Oz 


O(wt,,)  O(wt,,) _ O(wt,,) 
+ + + 


Pe yr age ne oe 


The two-dimensional compressible viscous flow energy equation 


reduces to: 


p— d SZ le+Sv*) +2 (pu) +e (pv) 
(27) 


Using the continuity equation and setting 
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E= (e+iv’)p, (28) 


the total energy per unit volume, the energy equation becomes 


SE 2 ((e+pul +S [esp vI 
: s (29) 
eae, +VT ~4.) + (ur +vt_-q,) 
ae eae, xy ~ Ux Oy xy yy dy 
where Tye Taye and Ty, are as previously derived (see Eq. 19). 


D. CONSERVATION LAW FORM OF THE GOVERNING EQUATIONS 
The conservation law form of the governing equations for 
the two-dimensional viscous compressible flows can be written 


as [Ref. 14] 


60, OF , OG _ OR, AS 56) 


ot dx dy Ox oy 
where Q is the vector of dependent variables, F and G are the 
inviscid fluxes along the x and y directions, respectively, 


and R and S are the viscous fluxes along the x and y direc- 


tions, respectively. These terms are given by: 
p pu pv 
pu u?+ puv 
= Fe P P Gr 
pv puv pv?+p 
E (E+p)u (E+p) v 
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T xx Ty 
R= Sex 
ae Rvs 
UT + VT xy” Fx ut xy? VT yy Gy 


Non-dimensionalizing the equations with 





» x «  U * . 
x= uj=— p*=- —P— p*= f. 
- —.Va pV 
mee +. pp -__@ ° t 
FOS, ees = = t’= 
sa ‘ H. = ve LEV, 


where EL is the reference length, and setting 


pP.VL 
H. 


Re = (31) 


the Reynolds number based on the reference length, obtain the 


non-dimensional form of the conservation law: 


00, OF, 0G _ i , OR, as 


—— —— — = — (+=) (32) 


ot ox dy Re ox dy 


with non-dimensional variable Qand the inviscid flum terms (F 


and G) as before, and the viscous terms given by: 


Re c (33) 


ar 
(y-1)M?Pr Ox 
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S= i. (34) 


oT 
xy VT yy (y-1)M*Pr Oy. 


Here Pr = C, p/K is the Prandtl number, and M = {V{/a,. In 
order to enable solution of the governing equation for 
arbitrary geometries, the governing equations are expressed 
for a curvilinear generalized coordinate system. The 
curvilinear coordinate system (&,n) is linked to the Cartesian 


coordinate system (x,y) by 


beh x VA EY 4 n=n(x,y,t) (35) 
The curvilinear coordinate space (—,n) is referred to as the 
computational domain and is linked to the physical domain 
(x,y) by the non-zero Jacobian of the coordinates transforma- 


tion 


O(b,n) ~ 08 On _ OF On . 1 
0(x,y) “Ox rig “Oy = Ox Oy _ Ox oy (36) 
oF On on 08 


It can be shown that the governing equations [Ref. 18] retain 
the strong conservation law form for a generalized coordinate 
system. The two-dimensional strong conservation law form for 


a generalized coordinate system is: 


22 + #36 | Be SB, 38) (37) 


* 8 "Re ob 
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and 


where 


p 
= t|pu 38 
” 5 fy iat 
E 


pu pV 
1 1 
Fz ae G=— ‘On, (39,40) 
(E+p) u- pot 


and U and V are the contravariant velocity components along 


the — and n directions, respectively given by 


v= Br Bry S ve Meu Sls yg (41,42) 
The viscous terms are given by 
0 
Tax ge tay ge 
R25 Tay GET yy SE (43) 
[ (ut ,,+Vt 


mn OT, 0& 
xy” (y-1) M? Pr ‘Ox Ox 
+ (ut wt VT yy -_  #P of oT 08 


(y-1)m? Pr Oy oy) 
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on on 
x Gy * To Sy 
co on 


+T 
xY Ox yy Oy 


-__ bs OT) ON 
[(ut+vt5, (yli) MP Pr Oy) Ox 
OT) ON; 


+ (UT LL tVt 
Nay Vy (y-1)m?Pr Oy oy 
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(44) 





III. NUMERICAL APPROACH 


A. NUMERICAL PROCEDURE 

The strong conservation law form of the two-dimensional 
Continuity, Navier-Stokes, and Energy Equations in generalized 
coordinates presented in Chapter II provides a useable format 
for implementing a numerical solution technique. The numeri- 
cal method used for the integration of the governing equations 
is a finite difference numerical scheme based on the Beam- 
Warming algorithm [Ref. 15]. The viscous terms are retained 
in both directions in order to enable capturing intense 
viscous effects encountered in massively separated flow 
regions at high angles of attack. The turbulent stresses were 
modeled using the Baldwin-Lomax eddy viscosity model [Ref. 
16). An algebraic C-type grid (157x58) was used for the 
computations. The boundary conditions were treated explicit- 
ly, and the unsteadiness was imposed by the motion of the 


grid. 


B. THE BEAM-WARMING ALGORITHM 
The strong conservation law form of the two-dimensional 


compressible Navier-Stokes equations in the vector notation is 


as follows: 














dU, OE(U) , AF(U) , OV, (U,U,) | eV, (U, U,) 


Rg ENT ge EEN ge AEE 


ot Ox oy Ox ox (45) 
, BWA, Uy) | OW, (U, Uy) 
oy oy 


here the vector U, the non-linear inviscid terms E and F, and 


the viscous terms R and S are given by: 


p pu pu 
pu pur +p puv 
U= E(U) = F(U) = 
pu pu? +p 
E, (E, + p)u (E, + p)u 
10] 
$ u(2u,. zi Uy) 
R=V,tV2.= (46) 
u(u, +v,) 


u(y + vx) + 3 uu(Qu, —vy) + kT, 
0 
uytou 
S = W, +W, _ 2 x) 
ju(2Qvy —Uuy) 


wu(uy + vy) + 5 yv(Qvy — uy) + kT, 


The Beam-Warming numerica] algorithm is an implicit finite 
difference scheme where the solution is marched in time using 


the difference formula 


8,At 9 At @ 0, 
nr = n Cape a 
ANU * aa, de On te oe Oa ame 


(47) 
+0[(8,-+-8,) (At)24(At)?] 
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where A°u = u™!- u". By substituting Eq. (45) into Eq. (47), 


obtain 





any = 21 4t 2 ane + A"V, +A"V ye 2 A" n n 
Teor hoe 1 a) + 55 (CAE + A"W, + A"W,) 
At 3 0 
ea? [Re cerevrevns 2 cer ew sw] (48) 





0, ‘= ] 
+5 +0; A” 'U+0  (« -4-0,) (At)? + can| 


The delta terms are linearized using truncated Taylor series 


expansions, so that: 
n+ =, n OE n n+1 n 2 
E = £E + (a5) (u™1-Uu") +O[ (At) ?] (49) 


which can be rewritten as 


AE = [A] 2A"U+0[(At)?) (50) 


where [A] is the flux Jacobian matrix dE/dU given by 





| 
0 | =f i oO 7 8 
| l | 
cet 1- 
ue + Sut | (y—3)u I (y—})v ! (i—y) 
2 2 | | i 
[A] =— \ | | (51) 
HU —v | —u { 0 
yE,u iE l | | 
+ (1 = yuu? +?) | —Yt 41 Gu? + v2) (y — 1)uu —yu 
. [ae 2 | 
In a like manner, A" F can be linearized as 
A°F = [B] 7A%uU+0[ (At) 2] (52) 


where [8] is the Jacobian matrix O@F/dU given by 
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| | 
0 | 0 | —-1 0 
| | 
uv I —v | —u 0 
{ i | (53) 
at 3- ]1- 
(8) = Le ee 1 w- Du ! (y—3)u ly—y 
pare ee ! 
| 
yE,u | | E —1 
res —y)vo(u? +7) 5 (y—A)uu | ee a +u?) —y 
| | 


| p 2 


The viscous delta term A" V,(U,U,) is linearized by writing 


a n 
A"V, = (a) A"U + (ss) A"U, + O[(An) } 


(54) 
= [P]"A"U + [R]" A"U, + O[(AD7] 


= ((P] —(Re])" A" + 2 ([R]" A"U) + O[(AN] 


where [P] is the Jacobian eV ,/au, [R] is the Jacobian av ,/avu,, 


and [R, ] = @[R]/dx. These matrices can be written as 








The matrix for [P]-{R,] is obtained by assuming that p and k 


are locally independent of U. In a like manner, A’ W,(U,U) is 


linearized as 


Aw, = ({o] - [5,]) 242v+2( [5] 2A") +OL(At)*] (57) 


where 
| | | 
0 ; 0 | oO | 
Hy, : My | 0 | 0 
1 | [ | (58) 
[2] — [S,] =-—— -v($ ) [Po aoe ) | 
| | 
4 | | 4 
-# ($1) —u? py , Uy | » ($x) 0 
md | | yy 
and 
r | | 9 
0 r 0 | 0 | 
| i 
a a P= : 1 1 (59) 
] | 4 | 
[S]=>5 —z HU 0 | 3H - 
| ; 
y . k 1k 
A. EN at) 2B 1 Gul ($u-E)oie 
($u-4)s (. et eye ! Cy 3 ty) 4% 


The cross-derivative terms can be evaluated explicitly without 


loss of accuracy by noting that 
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A°y? = A™1y,+0[(At)?] 
(60) 


A®w® = A™4w,+0[ (At) 7) 


for a uniform time step At. By evaluating the cross-deriva- 


tive terms in this manner, the block tridiagonal form of the 


final equations is maintained. 
Substituting Eqs. (50), (52), (54), (57), and (59) into 


Eq. (48) yields 





jin +f 5 ap (4) - 17} + (Rd)"— 25 RY" 


+2 (1B) (0) + Is)" - 95 %isi| Anu 








At C) G) 61 
eters gto ne 
6, At -$ 
n-1 a n-1 n-1 
wit 1-5-8 :) ae 
where [I] is the identity matrix. In Eq. (61), expressions 
such as ; 
[2 ( {al - (P) + (R,1) 7} A270 (62) 
are equivalent to 
(63) 


(2 ( (a) - [P] + (R,)) 7420) 
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The left-hand side of Eq. (60) is factorized in the following 


manner: 











8,ACc a 2 _# n 
({I] + 130, (3, CLA] - (PD + (Rd) ys (R171) 
@,At, ee rene (64) 
x({I] + 178, [= ((B]- fel+{s,]) grace })A®U 


= LHS(Eq.60) + O[ (At) 3] 


and the final form of the Beam-Warming Algorithm becomes 
LHS (Eq. (64)) = RHS[Eqg.(61)] . (65) 
The partial derivatives in the algorithm are evaluated using 


second-order accurate central differences. 


The Beam-Warming algorithm is implemented in the following 





manner: 
Step 1: 
6,At a Fe) (66) 
C12) are. (= ((A} - (P}+(R,)) alors [R} °]) AU, =RHS(61] 
Step 2: 
(try + BS (2 ([B)- fo] +15,]) ts} )) A°U=-A%U, (87) 
1+6, oy y ey? = 
Step 3: 


ul = U%4A"U (68) 


In Step 1, A” U, represents the remaining terms on the left- 
hand side of Eq. (64). Equations (66) and (67) represent 


systems of equations which have a block tridiagonal structure. 


28 











For the two-dimensional compressible Navier-Stokes equations 
the block matrices have a dimension of 4x4. Central 
differencing is used for the second order space derivatives. 
Dissipation terms for numerical stability are added. 

Warming and Beam [{Ref. 15] have shown that the algorithm 
can be simplified by assuming that p is locally constant so 
that @p/dx=0, d@p/dy=0. Then [P]-[R,]=0 and (Q1-[S,]=0. 
Therefore, the algorithm for 6, = %, @,= 0, implying second 


order accuracy in time, obtains the following form: 


[2] +4 [lal *- (R12 APU, = RHS (69) 


(ri +42 (2 a} *-Zis,) "7 v= AY, (70) 


In the present work the viscous terms were treated explicitly 
in order to avoid the expensive computation of the matrices R, 
and S\: The Beam-Warming algorithm with explicit treatment of 
viscous terms in generalized coordinates (&,n) is written as 
At @an, @pn n+l 

C(t] + [44 +5,8 ]) AU = 


(71) 


At(-SF"-2G%3R425%) = RHS® 


The algorithm is further simplificd by approximately factoriz- 
ing the LHS(61) or LHS(71) operator in order to avoid integra- 
tion of the full two-dimensional operator. The factored form 


of the algorithm is 
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( [2] +4£8,a2,) AUP} = RHS(61) ” ‘ees 
72 


( (7) +458, By.) AUPE = AUT? 

In practice, implicit algorithms have stability limits due to 
nonlinearities. In addition, whenever discrete methods are 
used to compute high Reynolds number viscous behavior, small 
scales of motion appear which cannot be resolved by the 
numerics. These scales are brought about by the nonlinear 
interactions in the convective terms of the momentum equa- 
tions. In any finite discrete mesh the small scales which 
cannot be resolved, result eventually in inaccuracy and 
contamination of the long wavelength, large scale phenomena. 
In order to dampen the high frequency numerical effects caused 
by the poor resolution of the small scales, implicit and 
explicit numerical dissipation is added to complete the 
algorithm. The numerical dissipation terms introduce an error 
level that does not interfere with the accuracy and resolution 
of any physical effects. The dissipation terms used in the 
present work are the ones suggested by Ref. 15. The complete 
factorized form of the numerical algorithms with the implicit 
and explicit dissipation terms is 


[T+ (AE) (O,A"+ (Dims) g) 1x {T+ (AE) (8,3 2+ (Dino) 1 470 a 


= At(-0,F"-0,G7+0,R"+0,5 °-€,,,:D,) 
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C. THE BALDWIN-LOMAX EDDY-VISCOSITY MODEL 

The Baldwin-Lomax model [Ref. 16] is an algebraic eddy 
viscosity model for calculating two- and three-dimensional 
separated flows. As opposed to the classical boundary-layer 
approximation which assumes zero normal pressure gradient in 
the boundary layer, thereby neglecting the normal momentum 
equation, the Baldwin-Lomax Thin Layer model retains the 
momentum equations and makes no assumptions about the pressure 
gradient. The advantage of this model arises in application 
to high Reynolds number, separated turbulent flows, including 
reverse flow regions. The Baldwin-Lomax model is a two-layer 
algebraic eddy-viscosity model and avoids the difficulty of 
finding the edge of the boundary layer. The effects of 
turbulence are simulated in terms of an eaay viscosity 
coefficient y,, where p+p, replaces p (the physical viscosity 
coefficient) in the stress terms of the governing equations. 


In the model, ph, is given by 


(p.); ysy 
be = t’ inner Cc (74) 


Me) guter ey 


where y is the normal distance from the wall and y, is the 
smallest y where inner and outer values are equal. In the 


inner region, the Prandtl-Van Driest formula is used: 


(Be) inner * p1?|w| (75) 
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where 1 = ky[1 - exp(-y*/A*)] and |o| is the magnitude of the 


vorticity. In two dimensions, 
Gu dv 
w| = (ou 8v (76) 
[eo ay Ox 
and 





yr = Pwey . yPuty ¥ (77) 


By By 


The subscript w indicates wall values. 


For the outer region 


(u,) outer ™ KC os P Praveen (VY) (78) 


where K, Cas are constants and 


Fuaxe = Ymax F max for boundary layers, or 
Fuaxe = Cwk Ymay U? pyp/F way for wakes and separated boundary 
layers. 


The quantities yy, and Fy, are determined from 


Fly) = ylwl[1-exp(y*/A‘)] (79) 


The exponential term in this equation is set equal to zero for 


wakes. 
Fy is the maximum values of F(y), which occurs at a 
value of y = Ymy- The function Fuiegl Y) is the Klebanoff 


intermittency factor given by 


Fuso (¥) = [1+5.5 (Cte y 64-2 (80) 
Yuax 
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The quantity Up), is the difference between U = U|y.my and 


minimum total velocity at a fixed flow-wise station: 
Unyp = (vie YF Doo (fat joss (81) 


where the second term is zero except in wakes. 
In order to achieve agreement with Cebeci [{Ref. 17], the 


values determined for the constants are: 


A = 26 k = 0.4 
Co. = 1.6 K = 0.0168 
Cries = 0.3 Pr = 0.72 
Cy = 1.0 Pr, = 0.9 


D. BOUNDARY CONDITIONS 

The following boundary conditions were used in the 
numerical implementation. A non-slip non-penetration boundary 
condition in terms of the contravariant velocity components 
was used on the airfoil solid boundary. The unsteady pitching 
motion was accomplished by rotating the grid about the quarter 
chord point at pitch rates determined by the desired reduced 
frequency. The inflow boundary was placed approximately eight 
chord lengths away from the body surface. Therefore, 
freestream conditions were specified at the grid inflow 
boundary. Simple first order extr=polation was used for the 
flow variables at the outflow boundary. For the wake, 
averaging of the flow variables on the upper and lower 


surfaces of the C-grid was used. The velocity on the surface 
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is given by U=X=0z and V H2=-Ox. The contravariant velocity 
components for viscous flow solutions are set equal to 0O. 
Therefore, the non-slip condition in terms of the physical 


velocity components for moving grid is expressed by 
ake ve A a 
V. Jl-n, €,J/ 12 


E. GRID GENERATION 
1. Grid Generation Methods 
In order to effectively compute complex flowfields, the 
physical domain of interest must be discretized with a finite 
mesh. The requirements of an efficient numerical grid are 
[Ref. 18]: 


1. smooth grid lines so that the transformation derivatives 
(metrics) are continuous. 


2. grid point spacing which varies inversely with expecta- 
tion of large numerical errors. 


3. minimizing grid skewness to avoid large truncation 
errors. 


Several general grid generation techniques exist; among these 
the most common methods are: 
1. Complex Variable methods, where the transformations are 
at least partly analytic; this method is restricted to 
two-dimensions. 


2. Algebraic methods, usable in two- or three-dimensions. 


3. Differential Equation methods, usable in two- or three- 
dimensions. 
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The grid generation method used in this study utilized the 
algebraic grid generation technique because of the computa- 
tionai efficiency anl specd it provides. 

2. Algebraic Grid Generation Method 

The algebraic method used employs known, easily 
invertible, functions to map arbitrarily shaped regions (in 
this case, the airfoil contour) into a simpler computational 
domain. The airfoil surface is unwrapped to form a simple 
curve in the computational plane as in Figure 6. In the 
computational domain lines are first drawn normal to this 
curve, and the grid points in the normal direction are 
subsequently generated. In the case of airfoil grid genera- 
tion, consideration must be given for the clustering of grid 
points near the airfoil in order to adequately resolve the 
near surface viscous layers. An algebraic function is used to 
provide a uniform stretching normal to the airfoil in the 
computational plane. The resulting computational domain grid 
is wrapped back to the initial physical grid using inverse 
transformations. 

The grid generation program used in this study yielded 
a 157x58 C~type grid. The program is listed in Appendix C. 
Grids for the airfoils under consideration are displayed in 
Figures 7-9. Figure 7 displays the global grid including wake 
for the NACA 0012 airfoil. Figures 8 and 9 show the body- 
fitted grid in more detail for the two modified airfoils. 


Grid clustering near the airfoils is shown for resolving the 
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Figure 6. Airfoil Grid Unwrapping 


boundary layer flow in turbulent and separated flow condi- 


tions. 
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IV. RESULTS AND DISCUSSION 


A. INVESTIGATION METHOD 

Two-dimensional unsteady flows involving both harmonically 
oscillating and rapidly pitching (ramp pitch-up) airfoils were 
investigated. Studies for the oscillatory cases were conduct- 
ed using the NACA 0012 airfoil in order to enable comparisons 
with previously conducted experimental measurements. Flows 
over all three airfoils rapidly pitched from 0° to 30° angle 
of attack with the reduced frequencies and Mach numbers (shown 
in Table 2) were subsequently investigated. 


TABLE 2. RAPIDLY PITCHING AIRFOIL CASES, Re = 4 x 106 


AIRFOIL 0.3 Mach 0.4 Mach 
k=.01 k=.02 k=.01 k=.02 


NACA 0012 Xx 


NACA 0012-63 X 


NACA 0012-33 





Solutions for the harmonically oscillating NACA 0012 
airfoil were obtained for flows at a freestream Mach number of 
0.3, Reynolds number based on the root chord Re =4x10° and for 
two reduced frequencies of k=0.1 and k=0.2. The reduced 
frequency is defined as k=wc/2U, for the oscillatory case 
where a(t)=a, + a, sin(ot). In terms of nondimensional quanti- 


ties k is given by k=/2M,, and &(t)=a, wcos(wt) is the 
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instantaneous pitch-up rate which varies during the cycle with 
w=2kKM, for a unit chord length. The present computations were 
conducted for a variation of the angle of attack as a(t)=10°+ 
6sin(t). Experimental test conditions for the same freestream 
Mach and Reynolds numbers were for a(t)=10°+ 5sin(t). Because 
the airfoil stalled just at 15° for the experiment, the 
computed flow for the same conditions did not yield a 
hysteresis loop. Therefore, as McCroskey suggested, the 
oscillation was increased by one degree and a hysteresis loop 
was obtained. The computed lift behavior shown in Figures 10 
and 11 exhibits the well-known hysteresis loop of a harmoni- 
cally oscillating airfoil experiencing dynamic stall. The 
computation was initiated from a steady-state solution at a=5° 
and was carried out for two cycles. Figures 10 and 11 display 
the results of the second cycle from the two-cycle computa- 
tion. 

Flow solutions for a rapidly pitching airfoil were 
obtained by pitching the airfoil at a constant rate froma 
zero angle of attack and steady-state flow conditions to an 
angle of attack of 30° at the desired reduced frequency and 
freestream Mach nu: er. For the case where ramp motion was 
imposed, the reduced frequency k is given by k=&c/2U, where «& 
is the constant pitch-up rate. In terms of nondimensional 
quantities k=o/2M,, or w=2kM,=constant, and a(t)=a(a-a,) ot, 


where in this study a =0° and @=30° and w=&(t). A summary of 
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the computed results for the rapidly pitching airfoils is 
provided in Table 3. 


TABLE 3. PEAK LIFT COEFFICIENTS 
Deeg ep OO ee thee ee eel ee tee Tl 


| AIRFOIL Mach Reduced Peak Angle of | 
No. Frequency Cy Attack 
0.3 -O1 1.90 19.40° 
NOO12 .02 2.10 23.40° 
0.4 -O1 1.97 21.08° 
02 2.24 25.21° 
0.3 -O1 1.84 18.50° 
NO012-63 -02 2.09 23.00° 
| 0.4 .O1 1.93 20.60° 
| 02 2.22 25.00° 
0.3 -O1 1.65 15.80° 
| NO0012-33 -02 1.94 19.60° 
| 0.4 .O1 1.74 $7. 10° 
| .02 2.09 23.15° 





A general observation on the computed solution is that the 


modified NACA 0012-63 airfoil exhibited comparable lift 
behavior to the baseline NACA 0012 airfoil. A slight angle of 
attack versus lift curve shift was observed between the two 
airfoils at all reduced frequency/Mach number combinations. 
However, for the same flow parameters the NACA 0012-33 airfoil 
consistently underperformed the larger leading edge radius 
airfoils. The difference in the unsteady lift behavior with 
increasing angle of attack resulted from the different flow 
character at the leading edge region as will be shown later in 
detail. In the following sections, detailed comparison of the 
lift behavior of the three airfoils 'is presented at the 


various flow conditions. The flow characteristics and the 
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development and progression of the dynamic stall process are 
examined. The effects of reduced frequency and freestream 


Mach number are also discussed. 


B. LIFT BEHAVIOR 
1. Harmonically Oscillating Airfoil 

Although hysteresis loops characteristic .* an 
oscillating airfoil undergoing dynamic stall were observed for 
reduced frequencies of both 0.1 and 0.2, agreement with 
experimental results varied. At a reduced frequency of 0.1, 
the response agreed well with McCroskey's experimental results 
(Figure 10) during the upstroke. Maximum lift coefficient was 
1.52 at an angle of attack of 15.8°. The lift behavior 
continues to agree well with the experimental data during the 
upstroke of the hysteresis loop and during the initial part of 
the downstroke. As the downstroke continues and the flow 
reattaches, however, the numerical solution displays signifi- 
cantly greater oscillations in lift coefficient compared to 
experimental data, which were obtained as an average over 
several cycles of oscillation. 

At a reduced frequency of 0.2 (Figure 11), the numeri- 
cal solution exhibited a much smaller hysteresis loop than the 
respective experimental data. In this case, maximum lift 
coefficient in the numerical sclution occurred just prior to 


the downstroke, so that flow reattachment occurred ajJmost 
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immediately. Maximum lift coefficient at k=0.2 was only 
slightly higher than at the k=0.1 case. 

The lack of agreement with experimental data at a 
reduced frequency of 0.2 and during the flow reattachment 
process at k=0.1 may be indicative of the poor behavior of the 
eddy-viscosity model (which is suitable for steady flows). 
Higher values of the reduced frequency resulted in larger 
discrepancies from the measured lift values. 

2. Rapidly Pitching Airfoil 

Lift coefficient vs. angle of attack comparison of the 
three airfoils is presented in Figures 12-15 for the Mach 
number/reduced frequency combinations listed in Table 2. For 
the same flow parameters, all three airfoils have nearly the 
same lift curve slope until the onset of dynamic stall. Alony 
with the aforementioned angle of attack shift of the lift 
curves between the NACA 0012 and -63 airfoils after stall 
occurs, the peak lift coefficients for the NACA 0012 airfoil 
are slightly higher than for the -63 airfoil. The peak lift 
coefficient for the NACA 0012 occurred 0.2 to 0.9 degrees 
angle of attack higher than the peak C, for the ~63 airfoil 
for the same flow conditions. With both airfoils having the 
same leading edge radius, the small difference in performance 


may be attributed to the slightly thicker contouring of the 


forward part of the baseline NACA 0012 airfoil (Table 1). 
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The maximum lift coefficient obtained with the NACA 
0012-33 airfoil was lower than that obtained with the other 
two airfoils and occurred at 1.8 to 4.0 degrees lower angle of 
attack. At a reduced frequency of k=0.01 and freestream Mach 
number of 0.3, the initial leading edge vortex originated at 
an angle of attack of 12° for the -33 airfoil. Under the same 
flow conditions, formation of the vortex occurred at 17° and 
16.5° for the baseline 0012 and -63 airfoils, respectively. 
Similar differences occurred at all other flow conditions 
investigated. The smaller leading edge radius of the NACA 
0012-33 airfoil promoted earlier development of the dynamic 
stall vortex from the leading edge upper surface. 

As will be shown, the dynamic stall vortex originates 
as a result of tic combination of the accelerated flow over 
the leading edge with the boundary layer reverse flow behind 
the leading edge. An adverse pressure gradient is encountered 
by the flow as it passes the suction peak just downstream of 
the airfoil leading edge. This adverse pressure gradient 
causes flow deceleration just aft of the suction peak, leading 
eventually to the boundary layer separating at a critical 
value of the adverse pressure gradient. The momentum of the 
flow, however, is increased as the flow accelerates around the 
leading edge and opposes the tendency of the flow to separate. 
Eventually, depending on the airfoil shape characteristics, 
the boundary layer separates, reversed flow occurs, and 


combines with the freestream to form thé leading edge vortex. 
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As observed in this investigation, the size of the 
leading edge radius is of primary importance in determining 
the angle of attack at which the boundary layer separates and 
rolls to form the dynamic stall vortex. For the same flow 
parameters and angle of attack during pitch-up, the critical 
pressure gradient for flow separation occurred earlier on the 
airfoil with smaller leading edge radius (NACA 0012-33). 
Stated otherwise, at a given angle of attack during pitch-up, 
the adverse pressure gradient aft of the suction peak is 
greater for the smaller leading edge airfoil. 

The contouring of the forward part of the airfoil is 
of secondary importance in alleviating development of the 
adverse pressure gradient. The effect of contouring is 
indicated by the slightly earlier development of leading edge 
reverse flow on the NACA 0012-63 airfoil compared to the 
baseline 0012 airfoil. 

In summary, enlarging the leading edge radius and 
thickening the contouring of the forward part of the airfoil 
results in decreasing the adverse pressure gradient encoun- 
tered by the flow, with resulting delay in separation. This 
is shown graphically in Figures 16-20 where instantaneous 
streamlines are shown for the three airfoils at the same 17° 
angle of attack, 0.4 Mach, and k=0.01. Reversed boundary 
layer flow has just begun on the NACA 0012 airfoil aft of the 
suction peak, whereas initiation of the dynamic stall vortex 


has already occurred on the NACA 0012-63 airfoil which has 
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thinner contouring. The NACA 0012-33 airfoil with the small 
leading edge radius has a fully developed dynamic stall vortex 
and is 0.6° angle of attack prior to dynamic stall. A 


detailed discussion of the vortical flow development follows. 
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0012-63, M=0.4, k=0.01, a=17 

















Figure 19. 














Figure 20. Instantaneous Streamlines, NACA 0012-33, M=0.4, 
-O1, a=17° ; 
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C. VORTEX FLOW DEVELOPMENT 
1. Rapidly Pitching Airfoil 

The development of the vortical flowfield as the angle 
of attack is increasing follows a consistent, general sequence 
of events. This sequence of events requires examination of 
the variation of several flow quantities as angle of attack is 
increased. Figures 21-33 illustrate the flow development 
stages using the velocity vector field, the pressure, and 
vorticity field contours for the NACA 0012-63 airfoil at M=0.4 
and k=0.01. Figures 34-47 display the associated surface 
pressures and temperatures of the developing flow. 

Initially, smooth streamlined flow is observed over the 
airfoil leading edge (Figures 21-22). As angle of attack 
increases beyond a certain critical limit depending on 
pitching rate, freestream Mach number, airfoil shape, and 
Reynolds number, small reverse flow regions develop on the 
upper surface aft of the suction peak and forward of the 
trailing edge due to the adverse pressure gradients encoun- 
tered by the flow. Figure 23 shows the reversed houndary 
layer flow near the airfoil surface. Figure 24 displays the 
developing leading edge vortex. 

The developing vortex forms on the upper surface just 
aft of the leading edge as a result of the combination of the 
accelerated flow near the suction peak and the reverse flow 
i:et aft of the suction peak. It is observed that the dynamic 


stall vortex initiates as a separated flow bubble, and rapidly 
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grows in size to form the characteristic leading edge dynamic 
stall vortical structure as the angle of attack increases. 
During this initiation process the primary vortex has a oval 
shape and becomes more rounded as the angle of attack increas- 
es. Figure 25 shows the oval dynamic stall vortex taking 
shape. Figure 26 displays the relative strength of the 
developing vortex. 

The vortex grows in size and strengthens, becoming 
approximately circular as angle of attack increases, and moves 
away from the airfoil. Figure 27 shows the development of the 
vortex as it moves downstream. Passage of the dynamic stall 
vortex over the airfoil surface induces reverse velocities and 
significantly contributes to the development of reversed 
flows. A secondary vortex originates near the leading edge 
and also grows as angle of attack increases. Figures 28 and 
29 show the development of a secondary vortex as the primary 
vortex moves downstream and promotes reverse flow over the 
whole upper airfoil surface. 

At the trailing edge, flow separation is initiated 
at approximately the same time as at the leading edge, anda 
pressure gradient exists between the lower and upper surfaces. 
The upper surface reverse flow combines with the high pressure 
flow from the lower surface to form a counter-clockwise vortex 
at the trailing edge shortly before dynamic stall occurs. 
Figures 30 and 31 show the developed trailing edge vortex 0.6° 


angle of attack before peak lift coefficient is obtained. As 
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the primary vortex grows and moves downstream past approxi- 
mately 60-70% chord, dynamic stall occurs, usually accompanied 
by development of a small, tertiary vortex at the leading 
edge. Figures 32 and 33 display the position of the dynamic 
Stall vortex 0.4° angle of attack after attainment of peak 
lift coefficient. 

By this time, the primary vortex is centered at a 
distance greater than the airfoil maximum thickness from the 
airfoil and combines with lower surface high pressure flow to 
energize the counter-clockwise vortex at the trailing edge. 

With further increase in angle of attack, the airfoil 
enters the deep stall regime, and the freestream flow has no 
effect on the upper surface flow characteristics. In this 
case, the aerodynamic behavior of the airfoil is greatly 
determined by the vortical flow field formed by the shed 
vortices, and the strength of the vortices themselves deter- 
mine the flow over the upper surface. Figures 48 and 49 show 
the NACA 0012-33 airfoil in deep stall at 0.3 Mach, k=0.02, 
a=22.0° and demonstrate that the vortices can combine to form 
a large vortical region and entrain other secondary flows. In 
this case a counterclockwise vortex was formed between the 


paired primary and secondary vortices and the airfoil. 
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Figure 21. Velocity Field, NACA 0012-63, M=0.4, 
a=14.0° 








Figure 22. Vorticity Contours, NACA 0012-63, M=0.3, k=0.01, 
a=14.0° 
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Figure 24. Vorticity Contours, NACA 0012-63, M=0.4, k=0.01, 
a=16.0° ; 
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Velocity Field, NACA 0012-63, 


Figure 25. 
a=17.0° 
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Figure 26. Vorticity Contours, NACA 0012-63, M=0.3, k=0.01, 
a=17.0° 
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Velocity Field, NACA 0012-63, M=0.4. k=0.01, 


Figure 28. 
a=19.0° 











Figure 29. Vorticity Contours, NACA 0012-63, M=0.4, k=0.01, 
a=19.0° . 
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Figure 30. Velocity Field, NACA 0012-63, M=0.4, k=0.01, 
a=20.0° 











Figure 31. Vorticity Contours, NACA 0012-63, M=0.3, k=0.01, 
a=20.0° 
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Figure 32. Velocity Field, NACA 0012-63, M=0.4, k=0.01, 
a=21.0° 














Figure 33. Vorticity Contours, NACA 0012-63, M=0.4, k=0.01, 
a=21.0° ; 
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Figure 34. Surface Pressure, NACA 0012-63, M=0.4, k=0.01, 
a=14.0° 
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Figure 35. Surface Pressure, NACA 0012-63, M=0.4, k=0.01, 
a=16.0° . 
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Figure 36. Surface Pressure, NACA 0012-63, 
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Figure 37. Surface Pressure, NACA 0012-63, M=0.4, k=0.01, 
a=18.0° ; 
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Figure 38. Surface Pressure, NACA 0012-63, M=0.4, k=0.01, 
a=19.0° ; 
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Figure 39. Surface Pressure, NACA 0012-63, M=0.4, k=0.01, 
a=20.0° 
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Figure 40. Surface Pressure, NACA 0012-63, M=0.4, k=0.01, 
a@=21.0° 
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Figure 41. Surface Temperature, NACA 0012-63, M=0.4, k=0.01, 
a=14.0° 


82 


4*10**6, k=0.01, a=16.0 deg 


— 


oO 
faa 
vt 
© 


Chordwise Position, x/c 





N0012-63 Airfoil, M 


\s : ie ‘© = 
— a S = 


WPINSIDI 0} VANLIIY sinjersduray, sovjans 





Figure 42. Surface Temperature, NACA 0012-63, M=0.4, k=0.01, 
“a=16.0° 
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Figure 43. Surface Temperature, NACA 0012-63, M=0.4, k=0.01, 
a=17.0° , 
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Figure 44. Surface Temperature, NACA 0012-63, M=0.4, k=0.01, 
@=18.0° 
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Figure 45. Surface Temperature, NACA 0012-63, M=0.4, k=0.01, 
a=19.0° . 
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Figure 46. Surface Temperature, NACA 0012-63, M=0.4, k=0.01, 
a=20.0° 
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Figure 47. Surface Temperature, NACA 0012-63, M=0.4, k=0.01, | 
a=21.0° : 
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Velocity Field, NACA 0012-33, M=0.3, k=0.02, 
22.G* 





Figure 49. Vorticity Contours, NACA 0012-33, M=0.3, k=0.02, 
a=22.0° 








2. Harmonically Oscillating Airfoil 


The process of vortical flow development at the leading 
edge during the upstroke is similar to that of a rapidly 
pitching airfoil examined in the previous section. Although 
the pitching rate is not constant in the case of the harmoni- 
cally oscillating airfoil, boundary layer separation and 
vortex development and shedding also occur during the up- 
stroke. 

The flow reattachment process during the downstroke 
involves the shedding of the large primary vortex into the 
wake and the diminishing intensity of the trailing edge 
vortex. Prior to the downstroke, the combination of the large 
clockwise primary vortex and the counter-clockwise trailing 
edge vortex cause extensive reverse flow, even outside the 
boundary layer, over the airfoil upper surface aft of 30% 
chord (Figure 50). As the angle of attack decreases further, 
flow over the leading edge upper surface rapidly reattaches, 
while the primary vortex is centered above the trailing edge 
(Figure 51). With further reduction in angle of attack, the 
primary vortex moves downstream of the airfoil, and the 
trailing edge vortex diminishes (Figure 52). Reverse flow at 
this time exists only on aft portions of the airfoil. As the 
angle of attack approaches 10-11°, the primary vortex has been 
swept downstream, the trailing edge vonecn has diminished, and 


£ 
4 


smooth, attached flow is established over the upper surface 


(Figure 53). 
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Figure 50. 


92 











Figure 51. Velocity Field, NACA 0012, M=0.3, k=0.1, «=13.0° 
downstroke 
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Velocity Field, NACA 0012, M=0.3, k=0.1, a=11.0° 


Figure 52. 
downstroke 
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Figure 53. 
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Velocity Field, NACA 0012, 
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D. REDUCED FREQUENCY EFFECT 

Al a given Mach number the three airfoils exhibited higher 
peak lift coefficients at respectively higher angles of attack 
at a redvced frequency of k=0.092 when compared to tnat of 
k=0.01. Lift curve slopes for a given airfoil at the two 
reduced frequencies were nearly identical to the point of 
dnitia* flow separation and subsequent dynamic stall at 
<=0.01. Lift coefficient at k=0.02 continued increasing, 
showing the increase in unsteady lift with higher pitching 
rate. The effect of increasing the reduced frequency at a 
freestream Mach number cf 0.3 is shown in Figures 54 and 55 
for the NACA 0012 and NACA 0012-33 airfoils, respectively. At 
the higher reduced frequency oi 0.02, formation of the primary 
vortex occurred at angles of attack 0.7° to 1.6° higher than 
at a reduced frequency of 0.01 in all cases. As a result, the 
entire vortical flow development was delayed to progressively 
higher angles of attack. Dynamic stalling angle occurred 4.0° 
to 4.8° higher at the higher reduced frequency and ihe 
resulting peak lift coefficients were 15-20% higher compared 
fo the k=0.01 cases (Table 3). 

Experimental work by Chandrasekhara and Carr [Ref. 3] 
using the NACA 0012 airfoil demonstrated the same trends, 
wherein a higher redced frequency resulted in retaining the 
vortex above the airfoil to higher argles of attack, thereby 


delaying dynamic stall. 
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Figure 55. 
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The effect of the increased pitching rate in delaying the 
onset of dynamic stall and increasing the unsteady lift may 
prove to be of practical operational benefit as knowledge of 


the details of the flow development increases. 


E. EFFECT OF FREESTREAM MACH NUMBER 

The effect of increasing freestream Mach number from 0.3 
to 0.4 resulted in slightly displacing the lift curve upward 
for the higher Mach number. The slope of the lift curve was 
unchanged. However, for a given angle of attack, a very 
slight increase in lift coefficient was observed at the higher 
Mach number before the onset of dynamic stall. Figures 56 and 
57 show the effect of increasing Mach number for the NACA 0012 
and 0012-33 airfoils at a reduced frequency of k=0.02. Flow 


at' a Mach number of 0.4 resulted in slightly higher peak Jift 


coefficients than at M=0.3. In addition, the peak lift 
coefficient at M=0.4 occurred 1-2° higher angle of attack than 
at M=0.3. 


The initial development of reverse flow regions and 
vorticity was largely independent of Mach number for a given 
airfoil in that the primary vortex originated at approximately 
the same angle of attack for the two Mach numbers. Subsequent 
growth of the primary vortex and development of the secondary, 
tertiary, and trailing edge vortices are delayed to slightly 
higher angles of attack at M=0.4 compared to the M=0.3 case. 


As a result, dynamic stall occurs at a higher angle of attack. 
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Figure 56. 
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El:raterinaris [Ref. 11] found that for the SSC-a90 airfoil 
with 9% thickness ratio, increasing Mach number from 0.2 to 
0.4 resulted in decreasing both lift coefficient and dynamic 
stalling angle of attack at a reduced frequency of 0.01. The 
experimental work by Chandrasekhara and Carr [Ref. 3] using a 
harmonically oscillating NACA 0012 airfoil also showed that 
increasing Mach number reduces the angle of attack at which 
dynamic stall occurred. Further experimental and computation- 
al investigation is required at comparable pitching rates to 
ascertain the accuracy of the computational model used at 
compressible flow Mach numbers. Further investigation into 
the influence of airfoil thickness ratio on Mach number 


effects is needed. 


F. PRESSURE GRADIENT AT FLOW SEPARATION 

The peak pressure gradient observed at initial flow 
separation occurring aft of the leading edge was investigated 
for the three airfoils at the reduced frequencies and Mach 
numbers listed in Table 2. Figure 58 displays a plot of 
streamwise pressure gradient along the airfoil chord for the 
NACA 0012-63 airfoil at different freestream conditions. The 
figure shows that, for the same airfoil, the peak pressure 
gradient encountered by the flow at initial flow separation is 
independent of freestream speed or pitching rate. Similar 


results were obtained for the other airfoils. Figure 59 shows 


the streamwise pressure gradient for the NACA 0012-33 airfoil. 








Figure 60 displays pressure gradients at initial flow 


separation angle of attack for the three airfoils at the 
freestream condition of 0.4 Mach and reduced frequency of 
k=0.01. The peak pressure gradient observed at the instant of 
flow reversal was higher for the NACA 0012-33 airfoil than for 
the airfoils with larger leading edge radius. Similar results 
were obtained at other freestream conditions. The streamwise 
location of the peak pressure gradient for the NACA 0012-33 
airfoil was slightly upstream of that for the NACA 0012-63 
airfoil as shown in Figure 61. The higher peak pressure 
gradient observed for the NACA 0012-33 airfoil is an indica- 
tion that flow momentum near the surface is higher around the 
smaller leading edge radius. The higher flow momentum 
combined with the relatively upstream location of the peak 
pressure gradient and sharper flow turning angle combinec to 
cause the magnitude of the peak pressure gradient for the NACA 
0012-33 airfoil to be 33% greater than that of the NACA 
0012-63 airfoil. Thus, the critical pressure gradient 
required for flow reversal is dependent on leading edge 
radius, in that a smaller leading edge radius increases the 
peak pressure gradient at flow reversal. 

tne streamwise location of flow separation is, in tvrn, 
dependent on the location of the peak pressure gradient. The 
actual location of flow reversal on the airfoil surface 
occurred 1.0-1.5% chord length downstream of the location of 


the peak pressure gradient as shown in Figure 61. This delay 
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is attributed to the time lag in the aerodynamic response to 
the rapid increase in pressure gradient. The more downstream 
location of the peak pressure gradient of the NACA 0012-63 
airfoil caused a similar relative location of the position of 


initial flow reversal compared the NACA 0012-33 airfoil. 


104 


4*10**6 


M=0.4, k=0.01, Re 





-~- 20012-63, a=16.0 
-- 00012-33, a=11.0 


epee a 


2 
~ 
= 
J 
a 
a 
Cc 
o 
=) 
z 


own OO NY OD NN O&O 
eg SES GL SD 


JUIIPVIH BNSsaig IIMUVIIS PIZIVULION 


Figure 60. 


Pressure Gradient Comparison, M=0.4, k=0.01 
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A. 


B. 


Vv. CONCLUSIONS AND RECOMMENDATIONS 


CONCLUSIONS 


The results of the unsteady flow so’utions displayed good 
agreement with the experimental] re vlts for the harmoni- 
cally oscillating airfoil at a reduced frequency of 0.1 
and Mach number of 0.3. At a higher reduced frequency of 
0.2, however, *ne agreement with experimental results was 
poor during the downstroke. 


The leading edge geometry of an airfoil was found to have 
a significant effect on the development of the vortical 
flowfield and dynamic stall characteristics of the 
airfoil. Of primary importance is the size of the 
leading edge radius. A larger leading edge radius delays 
development of the adverse pressure gradient necessary 
for boundary layer separation and eventual vortex forma- 
tion. Of secondary importance is the contouring of the 
airfoil aft of the leading edge. Thicker contouring 
forward of *he location of maximum thickness contiibutes 
to delaying flow separation. 


The effect of increasing the pitching rate of the airfoil 
is to enhance unsteady lift by delaying vortex formation 
to a higher angle of attack. The flow solutions present- 
ed in this study on the effects of reduced frequency are 
in good agreement with trends observed in experimental 
work. 


The streamwise pressure gradient required for flow 
separation is a function of airfoil leading edge radius 
and is independent of reduced frequency or freestream 
speed. 


RECOMMENDATIONS 


Conduct further experimental and computational investiga- 
tions to compare Mach number effects on flow development 
and dynamic stall. The investigations should be conduct- 
ed at comparable Mach and Reynolds numbers. 


Investigate the effect of airfoil thickness ratio on 
dynamic stall. 
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The understanding of harnessing the unsteady lift 
generated by rapidly pitching airfoils is only part of 
the knowledge required for practical implementation. 
Further investigations must include examining pitching 
moment development and how to minimize the adverse 
pitching moments produced. 


Conduct further computational work at higher reduced 
frequencies and Mach number to check the validity of 
using this eddy~viscosity model at those flow conditions. 


~J 


10. 


LIST OF REFERENCES 


McCroskey, W. J3., The Phenomenon of Dynamic Stall, 
National Aeronautics and Space Administration Technical 
Memorandum 81264, March 1981. 


McCroskey, W. J., and Pucci, S. L., Viscous-Inviscid 
Interaction on Oscillating Airfoils in Subsonic Flow, 
American Institute of Aeronautics and Astronautics Paper 
81-0051R, January 1981. 


Chandrasekhara, M., and Carr, L., Flow Visualization 
Studies of the Mach Number Effects on the Pynamic Stall 
of an Oscillating Airfoil. American Institute of Aeronau- 
tics and Astronautics Paper 89-0023, January 1989. 


Chandrasekhara, M. and Brydges, B., "Amplitude Effects on 
Dynamic Stall of an Oscillating Airfoil,” paper presented 
at the Aerospace Sciences Meeting, 28th, Rens, NV, 9-12 
January 1990. 


Mehta, U. B., Oynamic Stall of Osciilating Airfoils, 
Advisory Group for Aerospace Research and Development 
Paper #23 CP-227, Unsteady Aerodynamics, September 1977. 


Wu, J. C. and others, "Dynamic Stall of Oscillating 
Airfoils," paper presented at 42nd Annual Forum of 
American Helicopter Society, Washington, D.C., June 1986. 


Beddoes, T. S., Prediction Methods for Unsteady Separated 
Flows, Advisory Group for Aerospace Research and Develop- 
ment Paper #15 CP-679, 1980. 


Jang, H. M., Ekaterinaris, J. A., Platzer, M. F., and 
Cebeci, T., "Essential Ingredients for the Computation of 
Steady and Unsteady Blade Boundary Layers," ASME Paper 
90-GT-160, June 1990. 


Sankar, N. L., and Tang, W., Numerical Soluticn of 
Unsteady Viscous Flow for Rotor Sections, American 
Institute of Aeronautics and Astronautics Paper 85-0129, 
January 1985. 


Visbal, M. R., Effect of Compressibility on Dynamic Stall 
of a Fitching Airfoil, American Institute of Aeronautics 
and Astronautics Paper 88-0132, January 1988. 


1 


ml 


a 
dD 


1. 


Bes 


Ekaterinaris, J. A. Compressible Studies on Dynemic 
Stall, American Institute of Aeronautics and Astronautics 
Paper 89-0024, January 1989. 


Abbott, I. H. and Von Doenhoff, A. E., Theory of Wing 
Sections, Ch. 6, Dover Publications, Inc., 1959. 


Schlicting, H., Boundary Layer Theory, McGraw-Hill B3aok 
Company, 1979. 


Vivand, H., ‘Conservative Forms of Gas Dynamic Equa~ 
tions," La Recherche Aerospatiale, No. 1 pp. 65-68, 
January-February 1974. 


Beam, R. M. and Warming, R. F., "An Implicit Factored 
Scheme for the Compressible Navier-Stokes Equations," 
American Institute cf Aeronautics and Astronautics 
Journal, v. 16, no. 4, pp. 393-402, April 1978. 


Baldwi:.. B. S$. and Lomax, H., Thin Layer Approximation 
and Algebraic Model for Separated Turbulent Flows, 
American Institute of Aeronautics and Astronautics Caper 
78-257, January 1978. 


Cebeci, Teej Calculation of Compressible Turbulent 
Boundary Layers with Heat and Mass Transfer, Amer.can 
Institute of Aeronautics and Astronautics Paper 7O-7451, 
June 1970. 


Anderson, D. A., Tannehill, J, C., and Pletcher, R. H., 


Computational Fluid Mechanics and Heat Transfer, p. 490, 
Hemisphere Publishing Company, 1984. 


112 








APPENDIX A - USING THE PROGRAM 


A. THE NAVIER-STOKES CODE FOR GENERATING FLOW SOLUTIONS 

The main program NSE2D (listed in Appendix B for refer- 
ence) reads input data as supplied by the input file WNSE.1d, 
the grid files fort.11 or fort.21 as appropriate, and the flow 
file fort.31. The program forms the metrics and the Jacobian 
hy calling the METRIC subroutine. After calling the £LPS 
subroutine toa perform the computations, the main program 
outputs the rotated grid and the flow solution to the file 
fort. 20, and outputs the data of lift, drag, moment and 
pressure coefficients, time, angle of attack, and rotatiunal 
frequency to the loads file fort.32. 

The subroutine SLPS is the primary subroutine, and it 
calls other subroutines to perform the steps in the ADI 


algorithm. The major subroutines that SEZPS calls are: 


Subroutine Function 

MATRIX1,MATRIX2 Form block tridiagonal matrices for the & 
and y directions 

AMAT1, AMAT2 compute coefficient matrices o@F/dE— and 
dG/dAn 

STRESS supplies the viscous terms 

DISSIP adds fourth-order dissipation terms to 
the right-hand side of the Navier-Stokes 
Equations 


EXFPBC applies boundary conditions 





RES1 computes residuals from inviscid part of 
equations 
EDDY applies Baldwin-Lomax mode] 


B. DIRECTIONS FOR EDITING INPUT FILES FOR THE PURPOSE OF 

GENERATING FLOW SOLUTIONS 

1. Generating a Grid 

Initially, a grid must be generated as described in 

Chapter III in order to compute the flow about the airfoil in 
question. For this, it is necessary to use the AIRFGR.F 
program. With the rectangular coordinates of the desired 
airfoil, edit the input file AIRFGR.IN using the "ex" or the 
"vi" editors. Ensure that enough points are defined in 
reyions of high curvature such as the leading edge regions. 
Then ensure that no fort,.21 file exists by renaming or 
deleting it. At this time, run the grid program by typing 

AIRFGR < AIRFGR.IN 
After the program has run, ensure that the grid is satisfacto- 
ry by running FLOT3D. It may be necessary to define more 
coordinate points to ensure leading edge radius definition, 
upper or lower surface curvature, or trailing edge closure. 
Accuracy of trailing edge closure will determine wake thick- 
ness. The AIRFGR.F program is listed in Appendix C for 
reference. 

<. Rotating a Grid 

The grid generated in part 1 is for a zero degree angle 


of attack. Often it is desired to compute a steady flow or to 
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start an unsteady flow at an AOA other than zero. In order to 
do this the ROTGR.F program must be used. It is necessary to 
edit the ROTGR.F program and set the desired rotated grid 
angle by changing the al value in line 10. After the program 
has been edited, then: 

1. Compile ROTGR.F by typing "cf77 rotgr.f". 


“n 


Compiling ROTGR.F resuits in an output file named a.out. 
Move a.out to rotgr by typing "mv a.out rotgr". 


3. Execute the program by typing "rotgr". 


4. The output file produced by rotgr is fort.1?. It 
contains the desired grid rotation. Copy fort.1° to 
fort.11 or fort.21 as necessary for future flow selu- 
tions. 

The commands within the quotation marks should be typed, net 


> uuotation marks themselves. Upper and lower case are not 
important. 
3. Steady-State Solutions 


Ensure that the grid is rotated to the desired angle 


ef attack of the flow solution (part 2). The rotated grid 


should be copied to the file fort.21. The fort.21 file must 


heave the grid rotated to the desired angle of attack tor the 
steady-state solution! 
Modify the NSE.IN file as follows: 


line 2: 
~ set desired computational time interval (DT) 
- set ALFA, ALFAI, ALFA1 to zero 
- set REDFRE to desired reduced frequency 
- set AMINF to desired Mach number 


line 8: 
- set number of time steps (ISTP); usually 3000 are neces- 
sary for convergence at AOA#0O. In. order to run interac- 
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tively and not exceed time limit, set ISTP=1000 and run 
three times. 


line 10: 
- ensure that XREF=0.25 
- set TSHIFT=0.0 
- set REYREF to desired Reynolds number 
line 14: 
- For initial run, set RESTART OSCIL RAMP as FALSE TRUE 
TRUE. To restart (2nd and 3rd runs), set TRUE TRUE FALSE 
line 22: 
- Set time to 0.0 for 1st run. Set to preceding run time for 
restart. 
To run interactively and observe output as it develops, type: 
NSE < NSE.IN 
To run in the background, type: 
NSE < NSE.IN > NSE.OUT 
When the run is finished the last time step solution will 
consist of grid and flow solutions in file fort.20. It will 
be necessary to separate the grid and flow solution to restart 
the second run. The command "SPLIT20" will move the grid 
solution to fort.11 and the flow solution to fort.31, and will 
display the final run non-dimensional time. When restarting 
the second run, modify NSE.IN line 14 to RESTART=TRUE. A 
restarted solution uses grid file fort.j11 and flow file 
fort.31. Set time in NSE.IN to time specified for the 
previous run when SPLIT20 was executed. Execute the next run 
by the desired method as before. 


After convergence, save the grid and flow files 


(SPLIT20 to fort.11 and fort.31). 
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4. Unsteady Solutions 
Both oscillatory and ramp (constant pitch rate) 
solutions can be obtained by the WNSE.F program. One must 
ensure that the desired steady state initial angle of attack 
grid and flow solutions are obtained. The input file NSE.Jiv is 
modified as follows: 
line 2: 
- set time interval, freestream Mach number and reduced 


frequency, where k=&c/2U, for ramp solution and 
k=oc/2U, for oscillatory solution. 


- For ramp response: set ALFAI to steady-state i:ijitial 
conditions; ALFA, ALFA] are unused. 
For oscillatory response of form a=aa sin(ot), set 


ALFA=a-, ALFAI=a., ALFAI=(a--«:) since response will start 
at min e«. 


Jine 8: 
set ISTP 6000-12000 depending on length of response 
desired. 

line 10: 


set REYREF to desired Reynolds number. 

- set XREF=0.25 
eset TSHIFT = -0.5 so that response starts at min @ which 
is 1/4 cycle (-n/2 time shift as a=a, + aw. Sin(ot-n,?)} } 
from mean «a. : 


line 14: 
- set RESTART to TRUE 
~ set OSCIL to TRUE for oscillatory response (else FALSE) 
- set RAMP to TRUE for ramp response (else FALSE) 


line 19: 

- Program flow solution outputs can be selected for various 
angles of attack desired. Select desired AOA solutions 
by listing them multiplied by 100 (i.e. 11 degrees as 
1100). The outputs are listed as files fort.61 through 
fort.72. 

line <z. 


- Set time to 0.0 for initial run. For restart set time to 
the fort.20 output from the last timestep of the last run. 





Run the program by submitting it as a job to the queue 
by typing "qsub -1t 7200 subunst". . The subunst file is a 
command file which will order the steps for the Cray when the 
job number comes up in the queve. With at least 1 1/2 hours 
in computation, go run 5 miles, eat lunch, or play with the 
kids. Upon returning, check on queue status by typing "mqs”. 
If no entries, the run is finished and you can split the 
output files and continue, if desired. 

Any of the output files fort.61 through fort.7? and 
fFort.20 can be split into grid solutions and flow solutions 
for further analysis. Fort.20 file must be split to continue 
the response at further angles of attack, and to use the run 
time as the next run NSE.IN input in line 22. For continuous 
ramp or oscillatory responses, it is not necessary to reset 
ALFA, ALFAI, ALFA1, or RESTART; it is only necessery to reset 
desired AOA solutions to be recorded (line 19) and the time. 

5. Example 

As an example, to run an oscillatory solution for 
@=i04+47sin(t), it is necessary to get a steady flow solution at 
a=3° (3=10-7). The grid developed for the airfoil must be 
rotated to 3° using the ROTGR.F program and input to fort. 21. 
At this time a steady-state solution is obtained by modifying 


NSE.IN with the following inputs: 
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ALFA 0.0 

ALFAI 0.0 

ALFA1 0.0 

REDFRE as desired 

AMINE as desired 

TIME 0.0 

RESTART OSCIL RAMP FALSE TRUE FALSE 
XREF 0.25 

TSHIFT 0.0 


Run the program by typing " NSE < NSE.IN ( >NSE.OUT) " 

After 1000 timesteps, issue the command SPLIT20, so that the 
files fort.11 and fort.31 will be opened. Edit NSE.IN with 
the SPLIT20 time and set RESTART=TRUE. Run the program again, 
and then a third time, ensuring convergence (in lift coeffi- 
cient, drag coefficient, etc.) Save the final fort.11 and 
fort.31 files. 


Then, modify NSE.IN as follows 


ALFA 10.0 

ALFAI 3.0 

ALFA1 7.0 

ISTP 6000-12000 

TIME 0.0 

RESTART OSCIL RAMP TRUE TRUE FALSE 
XREF 0.25 

TSHIFT -0.5 

jal. ja2, ... ial2 as desired 


Run the program by issuing the command "qsub -lt 7200 
subunst". After completion, save the output files (fort.61 
through fort,.72) and fort.20. Do a SPLIT20, note time and 
AOA, then input TIME into NSE.IN for further runs. 

Graphs of the output files may be obtained by using the 
piot program PLOT3D. Observe output files (for example 


fort.65) and graph by typing the commands: "split" and then on 
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the 


and 


next line "65". Outputs will be in files fort.11 (grid) 
fort.31 (flow) which are read into the PLOT3L input. 


A typical NSE.IN file for the unsteady solution of the 


example proviem at k=0.02 and M=0.4 would be as follows. 


IMAX KMAX DT WW ALFA ALFA1L ALFAI REDFRE AMINE 
157 58 0.005 2.00 10.00 7.00 3.00 0.02 0.40 
ISPEC (FLAG FOR CHOOSING DIFFERENT SPECTRAL RADIUS) 
3 
WW2X, WW2Y, WW4X, WW4Y (EXPLICIT DISSIPATION COEF. FOR X AND Y) 
0.00 0.00 0.030 0.030 
ISTP NPER NOUT RES STRUNST 
9000. 18000. 1000. 100. 0. 
REYREF DMIN XREF TSHIFT 
4.00 0.00002 0.25 -0.5 
TSTARI 
-1.0 


RESTART,MULTIGRID OPTIONS SPECIFIED IN THE NEXT CARD 
TRUE TRUE FALSE 
CIRCOR ( CIRCULATION CORRECTION) 


TRUE 
31 127 

ITEL ITEU 

0300 0400 0500 0600 0700 0800 0900 1000 1100 1200 1300 1400 
ial ia2 ija3 ia4 ia5 ia6é ia7 ia8 ia9¥9 ialO iall ial2 
TIME 

0.0 
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SOURCE PROGRAM DATE 1/20/90 


nse2d.f TIME 5 1:46:04 am 


CHMSORRAKH SSAA E THERE AHH THOR E TE HESEREE AT SOURAL HHS OOREE THEORIES SE 
ce * 
ce MAIN PROGRAM . « 
qe « 
CUT Ae SSC He aR Perea eases ATHHHASAEEKHHSSHUAEEASOTARENHSSDHENKESSATEE HE 

PROGRAM NSFMAIN 

PARAMETER (1X=180,kx=60) 

COMMON /SURF/PSUR(IX) 

COMMON / FIX /OMEGA , BDOT 

COMMON /MUTUR/CMU (IX, KX) 

COMMON/SKINCF/CF(IX) 

COMMON /GR:D1/X(IX, KX}, Z(12%, KX) 

COMMON /PAR/GAMMA , REYREF , ALFA, ALFA1 , REDFRE , AMINF , ALFAI 

~OMMON/DGRID/!T.IMAX, KMAX, I TEL, ITEU 

COMMON /GRiD/YACOB( IX , KX} 

COMMON /DAMP/NW, WWI , WW2X, WH2Y ,WW4X, WHEY 

COMMON /FLOW/Q1 (14%, KX),Q2(1X,KXZ),Q3¢ IX, ER), Q4(IX, EX) 

COMMON/MTRIX/ XIR(IX, KX), X1Z( 1X, KA), ZETAR( 1X, KX), ZETAZ(IX, KX) 

L ,RITCIR, KX), ZETAT(IX, KX) 

COMMON /PLOT/TITLE (10) ,NSTPT, RESD( 3000) ,RES,CLE( 3000) ,CDPH( 3000) 

DIMENSION DRHO(IX, KX) 

COMMON/INITI/UINF, VINF 

COMMON /BC LOG/CIRCOR 

COMMON /LOGIC/RSTRT , PITCH, RAMP 

LOGICAL CIRCOR 

LOGICAL RSTRT, PITCH, RAMP 

CBARACTER ITITLE*8O 

COMMON/TITL/ITITLE 

COMMON/L2NORM/ RESDL2(10000) 

PI = 4. * ATAN(1.) 


a 


PROGRAM SOLVES TWO-DIMENSIONAL VISCOUS FLOW PAST ARBITRARY 
GEOMETRIES USING ADI PROCEDURE. 


° 
. 
e 
. 
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TAPES = FILE CONTAINING INPUT DATA 

TAPE6 = OUTPUT 

TAPES «= FILE TBAT SAVES THE FLOW FIELD AT THE END OF A RUN 
IP THE CURRENT RUN IS A RESTART OF A PREVIOUS RUN, THEN 
TAPE? IS LSED TO READ TRE FLOW FIELD INTO MEMORY 


READ INPUT DATA . 
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READ *, ITITLE 

READ (5,1) TITLE 

READ(5,20(1) 

READ (5,2) IMAX,RMAX,DT,WhK,ALFA,ALFA1,ALFAI , REDFRE, AMINF 


ISPEC « FLAG TO SPECIFY DIFFERENT SPECTRAL RADIUS FOR SCALING 
WITE TRE DISSIPATION 


READ(5,20€1) 
READ(5,2) ISPEC 
READ(5,20C1) 
READ(5,22i1) WW2X%,WH2Y, WHAX, HW4Y 
NSTP = NO. OF TIME STEPS TO BE DONE ON THIS RUN 
READ(5,20(1) 
READ (5,2:21) FNSTP,FNPER, FNOUT,RES 
NPER = FNFER 
NSTP = FNETP 
NOUT = FNCUT 


4 
4 
4 
4 
4 
5 
5 
5 
5 
5 
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REYREF= REYNOLDS NUMBER IN MILLIONS 

DNMIN = D:STANCE OF FIRST POINT OFF THE WALL 
FOR REYNOLDS NUMBERS UPTO 3 MILLION USE 0.00005 
MREF = RESERENCE VALUE AT X-AXIS 

READ (5,2001) 

READ (5,2221) REYREF,DNMIN,XREF,TSBIFT 

TSBIFT = TSHIFT * PL 

REREAL = "EYREF * 1000000. 

REYREF © REYREF * 1.£+06 


' 


SRO SII SENNA AAAAVAROD 


SPSTS 


A 
Na 
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TSTART = TIME THAT THE CALCULATIONS BAVE BEEN ADVANCED 
UPTO THE PFEVIOUS RUN. IF TSTART IS NEGATIVE THIS VALUE 1S 
OBTAINED FROM THE TAPE 3. 
READ (5,2001) 
READ (5,2221) TSTART 
2221 FORMAT(4F10.0) 
READ (5,2001) 
2001 FORMAT(1X) 
READ (5,2600) RSTRT,PITCE, RAMP 
READ(5,2001) 
READ(5,2000) CIRCOR 
FORMAT (3L5) 
read(5,*) ial ,ia2 ,ia3 ,ia4 , 1a5, 1a6, ia?, ia8, 1a9, ial10 
siail,ial2,1813,1a014,1a15 
falfed = 100.*alfad . 
IF( PITCH ) DT = PI / (NPER*REDFRE*AMINF) 
NEGATIVE REYREF MEANS INVISCID FLOW 


CONAvVawn—o 


C4 


PRINT OUT THE INPUT DATA 


ogee 
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WHI s 3. ° WH 
WRITE (6,4) TITLE 
WRITE (6,3) IMAX,KMAX,DT,WW,WWI, ALFA, ALFA1,ALFAI ,REDFRE, AMINF, XREF 
WRITE (6,€6) NPER, NSTP, NOUT 
FORMAT(/, 2%, SENPER= , 18, 5%, SENSTP= 18, 5%, SHNOUT=,18,/) 
IF(REYREF .GT.0.) WRITE (6,3700) REYREF 
WRITE(6,67) WW2X, WW2Y, WWAX, WHEY 
67 FORMAT(/,2X,5HNW2X=,F8.4,5X, SEWW2Y=,F8.4,52, SEWWAZ=,F8.4,5%, 
1 SEWW4Y*,F8.4,/) 
GAMMA=1.4 
ITEL = 31 
ITEU = 133 
TLE = ( ITEU - ITEL ) / 2 + ITEL 


CALL AZRFOL(ITEL, IPEU, ILE) 
IP(REYPEF .GT.0.) CALL CLUSTR(DNMIN)} 


J bj 
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READ(21) IMAX,KMAX 
read(21l) (( A(1.K),I"1,IMAX),R#1,KMAX ), 
(€ 2(1,R),T#1,IMAX),K=1,RMAX ) 


eda 


¢ 

CceeeC STARTING CONDITIONS. 

Cees DENSITY NORMALISED WITH RESPECT TO ROINE 

Cees VELOCITIES NORMALISED WITH RESPECT TO AINF 

Cere TOTAL ENERSY NORMALISED WITH RESPECT To (ROINP*AINF*AINE} 
c 


TOTEN=AMINF*AMINF*0. 5*+1./(GAMMA* (GAMMA-1.)) 
ALFA = ALFA * PI / 180. 
ALFAI * ALFAI * PI / 180 


Sireopeb-teraetet-tt- tet th drtetekoh $-act- tv rtd 
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SOURCE PROGRAM DATE 1/20/90 
nse2d.f 


ALMEAN = ALFA 
ALPAl = ALFAL * PI / 180. 
ALFAMAX = ALFA + ALMEAN 
ALFACR «© ALFAMAX 
ALFAS 15 THE ANGLE OF AIRFOIL WITH RESPECT TO FREESTREAM AT STEADY-STATE 
OR INITIAL POSITION OF UNSTEADY MOTION. 
IT SHOULD BE SET ACCORDING TO THE TYPE OF MOTION 
ALFAS = ALMEAK - ALFAL*COS(0.} 
UINF = AMINF * COS(ALFA) 
VINF © AMINF * SIN(ALFA) 
wnif = asinf 
winf = 0.0 
call rotgrid( imax,kmax,aifa ) 
DO 7 I=1,IMAX 
DO 7 Kel, KMAX 
Q1l¢1,K)"1. 
Q2(1,K)°UINF 
Q3(1,K)*VINF 
Q4(1,K)=TOTEN 
CONTINUE 
IF(RSTRT) THEN 
REWIND 11 
READ (11) IMAX , KMAX 
READ (11) ( ¢ X(I,K), Is1,IMAX }, Km1,KMAX 
( ( 2(I,R), T=l,IMAX ), K=1,KMAX 


gees 


a ee et 
rer 
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rewind 31 

READ (31) IMAX , KMAX 

READ (3]) FSMACH, ALFAD, REREAL, TIME 

READ (31) ( ( QLCI.K), Isl, IMAX ), K=1,KMAX 
( ¢€ Q2¢I,K), I=1,IMAX ), K=1,KMAX 
( . Q3¢1,K), Im1,IMAX ), K=1,KMAX 
« ¢ Q4¢1,K), Isl,IMAX ), Ke1,KMAX 


pel i a be, 


ENDIF 


if( time .gt. 500. ) time = 0. 

if( alfad .eq. 0 ) time = 0. 

TSTART = TIME 

IF(TSTART.GE.0.) TIME = TSTART 

IF(.NOT.(RSTRT)) TIME = 0. 

WRITE(6,62) TIME 

FORMAT(/ ,40BTHE CALCULATIONS STARTED AT TIME T* ,FL2.4,/) 
ASTART = ALFA + ALFA1 * SIN( 2*REDFRE*AMINF * TIME + TSHIFT ) 
ASTART = ASTART * ( 180. / PI ) 

IDONE = TIME / DT 

WRITE(6,6&1) ASTART, IDONE 

FORMAT( 2X. 1ORALFASTART=, F10.4,8X,15HITERATIONS DONE, 5X,16,/) 
CALL METRIC 

IMIN = ( ITEU - ITEL ) / 2 + ITEL 

ILOW = 2 * IMIN 

CHD = X(ITEL,1) - X(IMIN,2) 

NSTPT © NSTP + NSTPP 


We odatad aia ai ehas ae 4 
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‘3 


ome nna ——-> STARTING FROM A STEADY-SATE SOLUTION AT CERTAIN ANGLE <a> 
GIVE A TIME SHIFT <TSRIFT> TO START FROM Ao AND SET INITIALY 
TIME<0., KEEP TSHIFT THE SAME THROUGH THE GNST. CALCUL. 


aaan 


i 
i 
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IF ( PITCBR ) THEN 

DTP = PI / (NPER*REDFRE* AMINE) 
at = dtp 

endif 

alfacr = pis 9. 

alfamax= almean + alfal 

IF( RAMP ) READ (5,*°) TIME 

DO 1000 ITN=1,NSTP 

TIME = TIME + DT 


SANS 


REDUCE THE TIME STEP TO BALF AT THE PEAK OF THE CYCLE 


‘> Ob Gp & OD OP 


IF (PITCH) TBEN 


IF( ALFA .GT. ALFACR ) THEN 

DT = DTP * ( 1. ~ .5*(ALFA-ALFACR)/(ALFAMAX-ALFACR) ) 
ELSE 

DT = DIP 

END IF 
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OMEGA = 2.*REDFRE*AMINF * COS( REDFRE* 2.*TIME*AMINF +TSRIFT) 
1 *ALFA 

ALOLD * ALMEAN + ALFAL * SIN( 2. * REDFRE * AMINF * 
1(TIME ~ DT) + TSHIFT ) 

ALFA = ALMEAN + ALFAl * SIN( REDFRF * 2. * TIME * AMINF 
1 + TSHIFT ) 

ALFAD © ALFA * 45. / ATAN(1.) 

DALFA © ALFA - ALOLD 

CALL ROTGRID( IMAX , KMAX, DALFA) 

CALL METRIC 

END IF 

1F (RAMP) TEEN 

OMEGA » 2. © REDFRE * AMINF 

ALOLD © OMEGA * ( TIME - DT } 

ALFA = OMEGA * TIME 

ALFAD = ALFA * 45. / ATAN(1.) 

DALFA = ALFA - ALOLD 

CALL ROTGRID(IMAX, KMAX, DALFA) 

CALL METRIC 

END IF 

ALFAD © ALFA * 45. / ATAN(1.) 


rit pery 
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[eo oe OOOO 


DOSNACAaWN— OO SVNAVY EWN — OD RPNAV EWN —OOOVA 


CALL SLPS(ITN,ISPEC) 


APPLICATION OF BOUNDARY CONDITIONS 
CALL EXPBC(CL) 


CALL LOAD(CL,CDP,CDF,CM,ALFAS, XREF) 


Tyrtidd 
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PRINT OUT PRESSURE AT THE SURFACE 
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IF(ITN/50*S0.EQ.ITN) THEN 
WRITE (6,19) 
WRITE (6,33) ITN, TIME , DT 
IF(PITCH OR. RAMP) WRITE (6,3500) ALFAD, OMEGA 
IF{IIN/50° °500.80.2TN)CALL CPPLOT(ITEL, ITEU AMINE ,-X(1,1},2(1,1))} 
IF(ITW/407"400.8Q.1TN) MRITE(6,12) (1,CF(I),1=1,IMAX) 
WRITE (6,3000} CL , CDP , CDF , Ch 
END IF 
TIMPI = TIME / AMINE 
IF(1TN/50°50.£Q ITN)WRITE(3,8000)TIME, ALFAD,OMEGA,CL,CDP,CDF,CM 
IP(17W/1000°1000.£0.1T%) THEN 
DO 100 ° © IMIN , ITEU 


hid a 


an 
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SOURCE PROGRAM DATE 1/20/90 


nse2d.f 


eo 


BOC » { ACI,1) — R{IMIN,2) } / CBD 
WRITE(2,7000) 20C , PSUR(ILOW-1) , PSUR(T) 
MRYTE($, 7000) 20C , CE{ILOW-1) , CF(T) 
CONTINUE 

ENDIF 

CwU TO CALCULATE THE AVERAGE VALUES OF CL, CD, CM 

SCL = SCL + CL 

SCDF = SCDF + CDF 

SCDP = SCDP + CDP 

scm = SCM + CM 


DO 984 K = 1 , BMAX 
Do 986 J = 1, IMAK 
DRNO(I,K) = ABS{ Q1(1,K) - DRHO{I,K) ) 
IRES © IFIX(RES) 
IP(ITN/IRES*IRES EQ. ITN) THEN 
IPLOT « ( NSTPP + ITN ) / IRES 
RESD(IPLOT) = 0. 
CLE(IPLCT) = CL 
CDPE(IPLOT) = CDP 
po 985 R= 1, KMAX 
po 985 1 = 1 , TMAX 
RESD(IPLOT) = AMAX1(RESD(IPLOT) ,DRHO(I,K)} 
RESD(IPLOT) = RESDL2(ITN) 
ENDIF 
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IF(ITN/NOUT*NOUT .EQ. ITN ) THEN 
IF( ITN .£Q. 1*NOUT 10UT 31 
IF( ITN .EQ. 2*NOUT 1oUT 32 
IF( ITN .EQ. 3*NOLT IovuT 33 
IF( ITN .EQ. 4*NOUT IouT 34 
I¥( ITN .EQ. S*NOUT IOUT 35 
IF( ITN .EQ. 6*NOUT rout 36 
IF( ITN .EQ. ?*NoUT IovuT 37 
TB( ITN .FQ. B*NOUT IOUT 38 
IF( ITN .EQ. 9*NOUT IOUT 39 
IF( ITN .EQ. 10*NOUT 10UuT 40 
IF( ITN .EQ. 1i*NOUT rout 41 
IF( ITN .EQ. 12*NOLT 1ocT 42 
IF( ITN .EQ. 13*NOUT 10tT 43 
IF( ITN .EQ. 14*NOUT TOUT a4 
IF( ITN .EQ. 15*NOUT IOUT 45 
IF( ITN .EQ. 16*NOUT IOUT 46 


ann tune 


jout = 20 

REWIND IOUT 

WRITE (1OUT) IMAX , KMAX 

WRITE (IOUT) ( ( X(I,K), T=1,IMAX ), Re1,KMAX 
( ¢ ZCI,K), I=1,IMAX ), R=1,KMAX 

WRITE (XOUT) IMAX , KMAX 

WRITE (JOUT) AMINF, ALFAD, REREAL, TIME 

WRITE (IOUT) ( ( Q2CI,K), Ish, IMAX ), K=1,KMAX 
( ¢ Q2¢1,K), I=1,IMAX ), K=],RMAX 
( ¢ Q3¢1,K), I=1,IMAX ), Rel, RMAX 
( ( Q4¢1,K), I=l1,IMAX ), Kel, RMAX 

REWIND LOUT 


ENDIF 


dalfad * 100.*alfad 
fal = 100 
ia2 200 
ta3 300 
tad 400 
iad 500 
1a6 600 
ja? 700 
1a8 800 
fad = 900 
ialo= 1000 
iall= 1050 
dal2= 1100 
IF( IALFAD . . dal 
IALFAD . » is2 
IALFAD .EQ. ia3 
IALFAD . . iad 
IALFAC . . das 
IALFAD .EQ. 186 
IALFAD . 5 ia? 
TJALFAD . ia8 
IALFAD . » dag. 
IALFAD . » Lalo 
IALFAD . » dali. 
TALFAD . . ial2 
IF( IALFAD .EQ. ial 
IF( IALFAD . . a2 
IF( IALFAD .EQ. 1a3 
IE( IALFAD . . tea 
IF( IALFAD .EQ. ia5 
IF( IALFAD .£Q. 186 
IF( ITALFAD .EQ. 1a? 
IF( IALFAD .EQ. taé 
IF( IALFAD . . tags 
IF( IALFAD .EQ. ialo 
IF( IALFAD .EQ. 1ia11 
IF( IALFAD .EQ. 1a12 
REWIND IAOUT 
WRITE (ITAOUT) IMAX , KMAX 
WRITE (IAOUT) ¢ ( 2(I,R), Ie1,IMAX }, K=1,KMAX 
( ¢ 2¢C1.R), Yeh, IMAX ), K=1,KMAX 
(ITAOUT) IMAX , KMAX 
(TAOUT) AMINF, ALFAD, REREAL, TIME 
(IAOUT) ¢ ¢( QUCI.K), T=1,IMAX ). Kel, KMAX 
€ € Q2¢I,K), T#1,IMAX ), K=1,KMAX 
( ( Q3(I,K), T#l,1MAX ), K=1,RMAX 
( € Q4¢1,K), T=1,IMAX ), K=1,KMAX 
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REWIND IACUT 
ENDIF 
1000 CONTINUE 


REWIND 26 
(20) IMAX , RMAX 
(20) ( ( RCTLR), T#l,IMAX ), Kel, RMAX 
( ( 2C1,R), T#1,IMAX ), #1, RMAX 
(20) IMAX , KMAX 
(20) AMINE, ALFAD, REREAL, TIME 
(20) © ( QUCIL KR). T#l,IMAX ), R21, RMAX 
( Q2¢1.R), 1#1,TMAX ), R21, RMAX 
( O3CI.R), Tel, IMAX ), f=1,KMAX 
Q4(1.R), I=), IMAY ), RL, RMAX 








4 
12 
19 
22 
33 
2002 
3000 


3500 


‘3700 
: ¢5000 
6000 
7000 
8000 
8001 





SOURCE PROGRAM 





nse2d.f 


alfaft = alfa * (180/pi) 
write(6,667) alfaf, time 
format(17hANGLE OF ATTACK =,F15.16,5%,5HTIME=,F10.4,//) 


PRINT OUT VELOCITY PROFILE 


WRITE( 6,668) 

FORMAT(/', ,14X%, 2HUZ,142,28U2,10%, 8HDENSITY ,6X, 6BPRESSURE ,//) 
CMUL = AMINF / REYREF 

DO 4000 I = 70,130,2 

S= 0. 

po 4000 K = 2, 20 

S = S + SORT( (X(1,K)-Z(1,R-1))*92 + (21, R)-Z(1,R-1))**2 ) 
U = ( Q2(1,K)/Q1(1,K) ) 

Vi= ¢ Q3(I-R)/QL(I,E) ) 

UTOT = UfU + VeV 

P = (GAMMA-1.)*( O4(1,K) - .5*Q1(I,K)*UTOT ) 

WRITE (6,2002) I,K, U,V, Q1(I,K) , P 

CONTINUE 

CALL OUTRVC(REREFAL) 

STOP 

FORMAT(12A6} 

FORMAT (12A8 ) 

FORMAT(215,7F10.0) 

FORMAT (7F10.0) 

FORMAT (/, 5X, SHIMAX™ ,110,//, 5X, SEKMAX=,110, 

1/,5%.5B DT=,F20.8,/,5%,5HNW *,F20.8,/,5Z,5RNWI =,F20.8, 
2/,5X,5HALFA=,F20.8,/,5%,6BALFAL® ,F20.8,/, 5X, 6HALFAI= ,F20.8,/, 
+5%, 7EREDFRE= , F18.8,/, 5%, 6HAMINF™ ,F20.8,/, 5X, SEXREF= ,F21.8) 
FORMAT (1H1,5X,10A6) 

FORMAT (8(14,F10.4)) 

FORMAT (/) 

FORMAT (2F10.6,15) 

FORMAT(5), SRISTP=,15,5%, SBTIME=,F9.5,5%,3HDT=,F9.5) 

FORMAT (215, 5E14.6) 

FORMAT (5X, 3BCL=,F10.4,5%,48CDP=,F10.4,5X,4BCDP=,F10.4,5X,3HCM=, 
+F10.4) 

FORMAT (5X, SRALFA=, F10.4, 5%, 6HOMEGA* ,F10.4, 5X, 2HB=,F10.4, 5X, SHEDOT= 
1,F10.4) 

FORMAT ( 5X , 7HREYREF= , F20.4) 

FORMAT( 2110.6) 

FORMAT (8F14.8) 

FORMAT (3116.%) 

FORMAT (9113.6) 

FORMAT (33., BRCL(AVG )=,F12.7,4%, 9HCDF(AVG)*=,F12.7,4%, 9HCDP(AVG)=, 
+F12.7,4%,B8HCM(AVG)=,F12.7) 

END 
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SUBROUTINE AMAT1(K,IMX] ,XIX,X1Z,XIT) 
POCO OERATH OT Se RAH OO eH ERECT EAKESOETE TESST HOMERS CCE ENE RESeerREeAsete 
* 


c 
c 
ce SUBROUTINE AMAT1 © 
c 
¢ 


e . 
a 





ATSSSHSSHAKe Se SSERAKS PSHE RKASTSESHETKHO TSHR ARESESCHE TERRES eeREAEseecetenne 
PARAMETER (I1X=180,kx#60) 
\ COMMON/ FLON/O1 (IK, RX) ,Q2(1X, KX), Q3(1X, KX) ,04(1X, KX) 
! COMMON / PERTR/DQ1 (IX, KX) ,DQ2 (1X, KX) ,DQ3( IX, KX), DO4( 1X, KX) . 
COMMON/AM-A(4, 4,12) 
COMMON /PAR/GAMMA, REYREF, ALFA, ALFAL,REDFRE, AMINF , ALFAI 
DIMENSION 2IT(IX, KX) ,-XIX(IN, KX) ,AIZ(1IX, KX) 
REAL £0,£1,K2 


c 
Cees AMATI COMPUTES THE COEFFICIENT MATRIX DE/DQ DURING XI SWEEP 
c 
on = GAMMA - 1. 
{ Do 1000 1 = 2, ImMx1 
H Ko = XIT(1,K) 
: nl = XIX(I,K) 
; g2 = XIZ(1,K) 
d Uv = Q2(1,K) / O1(I,K) 
\ " = Q3(1,K) / Q1U(I,K) 
: EBSYR = O4(1,K) / QUCI.R) 
: PHI2 = 0.5 * GMl * (U * U+ We &) 
i TEETA = Fl *u+ 2° Ww 
: AQL,1,1) = KO 
: AQL,2.1) = Kl 
i A(1,3,.1) = K2 
438: A(1,4,1) = 0 
439! A(2,1,1) © Kl * PRI2 - U * THETA 
“aeo ! A(2,2,1) = KO + THETA - Kl © (GMI - 1.) * U 
44! A(2,3,3) = K2* U- GMl * Kl +8 
4425 A(2,4,1) = Rl * GMl 
442 A(3.2,13) = R2 0° PHI2 - W* THETA 
"444" A(3,2,T) = Kl * W- K2.* GMI * U 
445. A(3,3,1) = KO + THETA - K2 * (GM1 - 1.) * W 
446 | A(3,4,1) = K2 * GMI 
447, A(4,1,1) = TBETA * (2. * PHI2 - GAMMA * EBYR) 
448 A(4,2,1) = Kl * (GAMMA * EBYR - PHI2) - GML * U © THETA 
_ 449 | A(4,3,7) = K2 * (GAMMA * EBYR - PHI2) - GMl * W* THETA 
450. A(4,4,1) = KO + GAMMA * THETA 
451. 1000 CONTINUE 
452 RETURN 






















SOURCE PROGRAM DATE 1/20/90 
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mee Paproee 5k WEE 





SOURCE TEXT 
: SUBROUTINE AMAT2(I,KMU1,ZETAX, ZETAZ, ZETAT) 
455 | Coecccenne seer nse cceeennneeeeEnRseeeeNERACeCUNTKESSEH UNE HOE HERTS OEE 
456 | ce * 
457 ice SBBROUTINE AMAT2 £ 
458 ce * 
459  CoranseccrenserecesenreeceeneRneCRdeKEHOOOERRESSTOEURAAS TOS EEASSCSEEEED 
_ 460 | PARAMETER (IX=180,%x*60) 
461; COMMON /FLOW/Q1 (1X, KX) ,Q2(IX, KX) ,Q3(IX, KX), Q4( IX, KX) 
462 | COMMON/PERTR/DQ1 (IX, KX) ,DO2( IX, KX), DO3(IX, KX) ,DOA( 1X, KX) 
463 . COMMON /AM/A(4,4,T%) 
464 | COMMON /PAR/GAMMA,,REYRE F.ALFA,ALFAL,REDFRE, AMINE, ALFAI 
_ 465 | DIMENSIV. ZETAR( IY, 4X), 2ETAZ( IX, MQ), ZETAT( IX, KX) 
. 466 REAL KO,K1,%2 
467 Cc 
ys ag8 | cess AMAT2 COMUIES THE COEFFICIENT MATRIX DF/DQ DUEING ETA SWEEP 
469 5c 
470 | GMi = GAMMA - 1. 
4271) DO 1000 K © 2, KMXI 
472: 50 = ZETAT(I,K) 
473 Rl = ZETAX(1,K) 
474; x2 = ZETAZ(I,K) 
| 475} U = Q2(1,K) / QLU(I.K) 
476: w = Q3(I1,K) / QL(I,R) 
477 | £BYR = Q4(1,K) / QL(I.R) 
478, PHI2 = 0.5 * GMl * (U ®t + We w®) 
479 | THETA = KL *U+ k2* Wh 
480 | A(1,1.8) = KO 
48) | A(1,2,F) = Kl 
_ 482: A(1,3,5) = K2 
483 | A(1,4,K) = 0 
484, A(2,1,K, = Kl * PBI2 -U * THETA 
485 | A‘2,2,K) & FO * THETA - Kl * (GM1-1.) * U 
486 A(2,3,K) = K? * U ~- GMl * Kit W 
487 | A(2,4,R) = Kl * GM1 
488 A(3,1,K) = S2 * PHI2 - W + THETA 
489. A(3,2,K) 2 Kl * «© - K2* GML * U 
490. A(3.3,K) = KO + THETA -K2 * (GMI-1.) * B® 
491. A(3,4,R) = K2.*° GM] 
492 | A(4,1,8) = TRETA * (2) * PRI2 GAMYA * EBYR) 
493 AC4,2,8) © Kl * (GAMMA * EBYR - PEIZ) - GM] * U * TRETA 
494 A(4,3,K) = K2 * (GAMMA * EPYR - PHID) - OM1 + Wo* THETA 
495 - A(4.4,R) = KO + GAMMA * TRETA 
496 1000 CONTINUE 
497 RETURN 


498 | END 
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SOURCE TEXT 
SUBROUTINE SLPS(ITN, ISPEC) 


CEPR EAHHSS STE MAAS OHHH KHSHSHETEASHERERLEDESERRASSSSERAAHSOSERAESOEEE 
ce e 
ct SUBROUTINE SLPS . 
ce * 
CAAAAPHSHASTKH FS STAC AHSSHHARKSSSARRESSSORRAKSSSEERTATOTOLEKESESSERAED 

PARAMETER (IX=180,Rx=60) 

COMMON/FLOW/Q1 (IX, KX), Q2 (IX, KX) ,O3{IX, KA), Q4(1X, FX) 

COMMON/FIX/OMEGA , DOT 

COMMON/PERTR/DQ1 (1X, RX) ,DQ2( IX, KX) ,DQ3(1X, KX) ,DQA(IX, KX) 

COMMON/AM. (4,4,1%) 

COMMON/TRID1/DD(4,4,1%, KX} 

COMMON/TR2D2/MM(4,4,1X, KX) 

COMMON/TRID2/EE(4,4,1%, KX) 

COMMON/TRID4/GG (4, 1X, KX} 

COMMON/PAR/GAMMA , REYREF , ALFA, ALFA1 ,REDFRE, AMINF , ALFAL 

COMMON/DGRID/DT, IMAX, KMAX,ITEL, ITEU 

COMMON/GRID/YACOB( IX, KX) 

COMMON/DAMP/}W , WHI, WH2X , WN2Y , WWAX, WAY 

COMMON/MTRIX/ XIX(IX,EX),-IZ(IX, KX), ZETAX( IX, KX) , ZETAL(IX, KX) 

1 ,RIT(IX, RX), ZETAT( IX, KX) 

COMMON /RADIUS/ SPECX(IX, KX), SPECY (IX, KX) 

COMMON /L2NORM/ RESDL2(10000) 

REAL MM, DELTAT(IX, RX} 


Tul 


a 


it 
v 
WN =—QCOONAVARWN—ODONAVAWN =O 


{2 


edn eu 


ee ee eel ie ter eee ence Cnceb alee oy 


TRE SUBROUTINE SLPS DOES TYE BULK OF THE WORK FOR THE ADI ALGORITHM. 
IT CALLS FLUX AND COMPUTES YIGHT BAND SIDE DURING TRE TNO SWEEPS, 
ASSEMBLES THE COEFFICIENT MARICES, ADDS IMPLICIT AND EXPLICIT 
DISSIPATION AND CALLS THE TRIDIAGONAL SOLVER TO OBTAIN TRE FINAL 
SOLUTION. 

IMl = IMAX 

IM2 «© IMAX 
KM1l = KMAX 
KM2 = KMAX 

TBE DISSIPATION TERMS ARE CONSTRUCTED AND STORED IN THE ARRAYS 001, 
DQ2,DQ3 AND DOA, 


CALL SPECF (ISPEC) 
CALL DISSIP 


ON TO 0Q1,DQ2,DQ3 AND DOQ4. 
THE RIGHT HAND SIDE AT KNOWN TIME LEVEL IS NOW COMPUTED AND ADDED 
CALL RESI 


DO 8999 K 
bu 8999 I 

DELTAT(1,K ; 
IF((REDFRL.LT.0.001).AND.ITN.EQ.1) THEN 

DO 1K = 2, KML 

polis 2, i) 

DELTAT(I,R} = 1.0 / ( 1. + SQRT(ABS(YACOB(I,K)) ) ) 
ENDIF 


) 


IF VISCOUS FLOW IS COMPUTED THE VISCOUS TERMS ARE ADDED TO DQ] ETC. HERE. 


TF(REYREF./ 0.) CALL STRESS(ITN) 

I-SWEEP. 

DTH = DT + 0.5 

DTK = DT * Wal 

DO 3K = 2, KM 

CALL AMATI(K,IMAX-1,X1X,XIZ,XIT) 

po 411 14 

DO 4 12 

pe 5 I 

£E(11,12,3-1,%) A(11,12,141) * DTH * DELTAT(I,R) 
DD(11,12,1-1,K) A(IL,I2,I-1) * DTH * DELTAT(1,X) 
CONTINUE 

CONTINUE 


IMPLICIT DAMPING ADDED BERE. 


= 2, IM 
= SPECX(1,K) * DTW * DELTAT(I,K) 
© 2. * SPECX(I,K) * YACOB(I,K) * DTW*DELTAT(I,K) 
DD(1,1,I-1,K) - DTl * YACOB(I-1,K) 
DD(2,2,1-1,K) - DTl * YACOB(I-1,R) 
DD(3,3,1~-1,K) DT1 * YACOB(I-1,K) 
DD(4,4,1-1,R) DT] YACOB(I-1,K) 
£E.(1,1,1-1,K) DT1 YACOB(1+1,K) 
EE(2,2,1-1,K) DTl YACOB(I+1,K) 
££(3,3,1-1,%) boT1 YACOB(I+1,K) 
£E(4,4,1-1,K) pTl YACOB(1+1,K) 

+ + DTWI 

. * DTWI 

. * DTK 

. * DTWI 


CONTINUE 
CONTINUE 

2, KM 

2, 1M) 

DOQ1(I,K) * DELTAT(I,R&) 
DOQ2(1,K) * DELTAT(1,R) 
DQ}(1,K) * DELTAT(I,R) 
DQ4(1,K} * DELTAT(I,K} 


7eenne 


CONTINUE 
PERFORM BLOCK TRIDIAGONAL MATRIX INVERSION FOR TEE ENTIRE PLANE 


CALL MATRX) (IMAX, KMAX) 

DO 991 K 2, KM 

DO 991 TI 2, Tel 
DQ1(1,K) GG(1,1-1,R) 
DO2(1,K) GG(2,1-1,K) 
DO3(1,R) GG(3,I-1,K) 
DQ4(1,F) GG(4,1-1,K) 
CONTINUE 


K-SWEEP BEGINS HERE. 


po 131 * 2, IM 

CALL AMAT2(%,KMAX-1,ZETAX, ZETAZ. ZETAT) 

po 15 Ti 1,4 
1 
2 


bo 15 12 
po 15 kK 
EE(I1,12,1,8-1} A(T1,12,KR+1)*DTB * DELTAT(1,#} 
DD(11,12,1,R-1) -A(I1,12,K-1)*DTB * DELTAT(I,&) 
CONTINUE 


PARWMAAAQAGSG 
auavawn—o 


SECOND OR"NER DAMPING ADDED BERE. 











SOURCE PROCRAM 


nse2d.f 


DO 16 K = 2, EML 

DT1 = SPECY(I,K) * DTW * DELTAT(I,K) 

DTwI = 2. * sPECY(I YACOB(I,K)} * DTW*DELTAT(I,K) 
DD(1,1,1,K-1) DD(1,1,1,K-1) YACOB(I,KF-1) 

DD(2,2,%,R-1) = DD(2,2,1,R-1) YACOB(1,K-1) 
bD(3,3,1,K-1) DD(3,3,1,K-1) YACOB(I, K-12) 

DD(4,4,1,K-1) = DD(4,4,1,R-1) YACOB(1,K-1) 

BE(1,1,1,-1) £EE(1,1,1,K-1) TACOB(I,K+1) 

EE(2,2,1,R-1) EE(2,2,1,K-1) YACOB(I, K+1) 

£E(3,3,1,K-1) ££(3,3,1,8-1) YACOB(I, K+1) 

BE(4,4,1,f-1) EE(4,4,1,K-1)} YACOB(I,K+1) 
MM(1,1,1,K~1) + DTWI 

MM(2,2,1,K-1) 
hM(3,3,1,K-1) 
MM(4,4,1,R-1) 
CONTINUE 

po 17K 

bo 171 
GG(1,1,k-1) 
©G6(2,1,K-1) 
GG(3,1,K-1} 
GG6(4,1,K-1) 
CONTINUE 


+ DTWI 
. + DTWI 
+ DTwWI 


Pere 


2, EM. 
2,0 
DQ1(1,K) 
DQ2(1,K) 
DO3(1,K) 
DQ4(1,K) 


PERFORM BLOCK TRIDIAGONAL MATRIX INVERSION FOR THE ENTIRE PLANE 


CALL MATRX2(IMAX, KMAZ) 

DO 18 K =2, ML 

po 18 I 2, IM 
DQ1(1,R) GG(1,1,K~1) 
DQ2(1,K) GG(2,1,K-1) 
DQ3(1,K) GG(3,1,K~1) 
BQ4(1,KR) GG(4,1,R~1) 
CONTINUE 


UPDATE FLOW VARIABLES AT INTERIOR POINTS. 
CONTINUE 
RMAX 
RUMAX 
RVMAX 
EMAX 

DO 995 K 
po 19 1 
Ql(1,K) 
Q2(1,F) 
Q3(1,K) 
O4(1,F) 
CONTINUE 
DO 995 1 = 2, IM1 

IF (RMAX.LT.ABS(DQ1(1,R)*¥YACOB(1,K))) THEN 

IR =I 

KR = K 

END IF 

RMAX = AMAX](RMAX,ABS(DQ1({I,K) * YACOB(I,K) 
RUMAX = AMAX1(RUMAX,ABS(DQ2(1,K) * YACOB(1,K 
RVMAX = AMAX1(RVMAX,ABS(DQ3(1,K) * YACOB(1,K 
EMAX = AMAX1(EMAX,ABS(DQ4(1,K) * YACOB(1,K) 
CONTINUE 


DQ1(I,K) * YACOB(I,K) 
DQ2(1,K) * YACOB(I,K) 
DQ3(1,K) * YACOB(I,K) 
Q4(1,K) + DQ4(1,R) * YACOB(I,K) 


ev ene euengn 


COMPUTE L2 NORM OF Q1, Q2, Q3, AND Q4 RESIDUAL 


R12 = 0. 
RUL2 = 0 
RVL2 = 0. 
El2 = 0. 
Do 996 K= 2, KM1 
DO 996 IT * 2, IM1 

RL2 = R12 + DQL(1.K)**2 

RUL2 * RUL2 + DQ2Q(1,K)**2 

RVL2 = RVL2 + DQ3(1,K)**2 

El2 = El2 + DQ4(1,K)**2 

CONTINUE 

RL2 = RL2 / ( IM2*KM2 ) 

RUL2 = RUL2 / ( IM2*KM2 ) 

RVL2 = RVL2 / ( IM2*KM2 ) 

E12 EL2 / ( IM2*KM2 ) 

R12 SCRT(RL2) 

RUL2 * SORT(RUL2) 

RVL2 = SORT(RVL2) 

EL2 = SURT(EL2) 

RESDL2(ITN) = RL2 

TF((ITN-1)7500*500.£Q. (ITN-1)) WRITE (6,3002) 

IF(ITN/S0*50.EQ.1TN) WRITE (6,3001) RMAX,RUMAX,RVMAX, EMAX,IR,KR 

IF(ITN/50*50.EQ.ITN) WRITE (6,6004) RL2 ,RUL2 ,RVL2 ,EL2 

RETURN 
3002 FORMAT(/,5X,'DRMAX',11X%,'DUMAX',112%, 'DVMAX',11%,'DEMAX' , 9x, 

1'IR',3%,'KR’) 
3001 FORMAT(4(£14.6,2X),215) 
6004 FORMAT( 4(E14.8,2X) ) 











SOURCE PROGRAM 


SUBROUTIN® MATRX1( IMAX, KMAX) 


POSCHARAARASASHAERAHSESTERE ARES SATUERARACSSSARRASOSSCETAEHAESHEEAAAHSOCERERTOCER 
2 e 
SUBROUTINE MATRIL ¢ . 
* 
SHAS SOSHERAEH SE SHE TAASSHERTATSSSHERKHKSSHSELRKASSHERARKRATASSEREKASSCE EAE 
PARAMETER (IX=180,kx=60) 
COMMON/TRID1/DD(4,4,I%, RX) 
COMMON /TRID2/MM(4,4,IX,KX) 
COMMON, TRID3/EE(4,4,1X,KX) 
COMNON/TRID4/ GG(4, IX, KX) 
COMMOK, DTESTAC4 4,20}, 28(4 6 TR O20) ,0(4, 5,12} 
REAL MM 
REAL 111,121,131, 141,122,132,142,133,143,144 
/LLT,L21,L31,L41 


THIS SUBROUTINE PERFORMS THE BLOCK TRIDIAGOWAL MATRIX IVERSION FOR 
A ENTIRE PLANE DURING THE Xi- SWEEP 


EM1 © KMAX - 1 

pOltll=1, 4 

pOo1lK«= 2, RM 

AI = 1. / MM(1,1,1,K) 
GG(11,1,K) = GG(11,1,K) * AI 
BH(I1,1,1,6) = EE(I1,1,2,K) 
BH(I1,2,1,K) = EE(I1,2,1,K) 
WBE(11,3,1,K) = EE(11,3,1,K) 
BH(T1,4,1,K) = EE(I1,4,1,8) 
CONTINUE 


DO 1000 1 = 

Do S i= 

DO 2 Par 

C(11,1,K) DD{I1,1,1,K) * GG(1,1-1,K) 
1 DD(11,2,1,R) * GG(2,1-1,K) 
2 DD(11,3,1,K) * GG(3,I-1,R) 
3 DD(I1,4,1,K) * GG(4,1-1,R) 
2 CONTINUE 

pDO5 12 

pos «K = 

A(I1,12,K)= MM(11,12,1,%) 


~ DD(IL,1,1,K) *-BH(1,12,I-1,R) 
~ DD(I1,2,1,R) * BH(2,12,1-1,R) 
~ DD(11,3,1,K) * HB(3,12,1-1,%) 
~ DD(TL,4,1,R) * HH(4,12,1-1,R) 
) 


1 
2 
3 
C11, 22-1, K)= EE(I1,12,1,K 
5 CONTINUE 
DOJ K= 2, RMI 
Lll = a(1,1,8) 
Lll = 1. / Li 
Ul2 = A(1,2,R) 
V13 = A(1,3,R) 
U14 = A(1,4,R) 
L221 = A(2,1,R) 
L31 = A(3,1,K) 
L41 © A(4,1,R) 
L22 = A(2,2,R) - L221 * U12 
121. = 1. / 122 
U23 = (A(2,3,K) - L21 * U13) * L2I 
U24 = (A(2.4,R) - L22 * U14) * Laz 
L32 = A(3,2,R) - L31 * v12 
L42 = A(4,2,K) - Lél * Ul2 
L33 = A(3,3,K) - L31 * U13 - 132 * U23 
L3I = 1. / 133 
U34 = (A(3,4,K) ~ L31 * Ul4 - 132 * U24) * L3T 
L43 = A(4,3,R) - LS1 * UL3 - L42 * U23 
144 = A(4,4,R) - LAl * U14 - 142 * U24 - L43 © 034 
“4r= 1. , 144 
C(1,1,K) = C(L,1,K) * Lit 
C(2,1,R) = (C(2,1,R) - L221 * €¢1,1,K)) * L2I 
€(3,1,K) = (C(3,1,) - 31 * C(1,1,R) 
1 - 132 * €(2,1,K)) * L31 
C(4,1,K) = (C(4,2,R) - L41 * C(1,1,K) — L42 * €(2,1,R) 
1 - bdo + €(3,1,K)) * LAT 


C(1,2,R) = €(1,2,K) * LUI 

C(2,2,8) = (C(2,2,) - L21 * C(1,2,Kj} * L2r 

€(3,2,K) = (C(3,2,K) - L31 * C(1,2,R) 

1 - 132 * €(2,2,K)) * L32 

C(4,2,K) = (C(4,2,R) - L4l © C(1,2,FR) - L42 * C(2,2,K) 
1 - L43 * €(3,2,K)) * LAE 
C(1,3,R) = C(1,3,K) * LIT 

C(2,3,K) = (C(2,3,K) - L22 * C(1,3,K)) * L2I 

c(3,3,K) (C(3,3,K) - 132 C(1,3,K) 

1 - L32 * €(2,3,K)) * 132 

C(4,3,8) = (C(4,3,K) - L4) * C(1,3,R) - L42 * C(2,3,K) 
1 - L43°* €(3,3,K)) * LAI 
C(1,4,8) = C(a,4,R) * LI 

€(2,4,K) = (C(2,4,K) - L21 * C(2,4,K)) * L2I 

C(3,4,R) = (C(3,4,8} - L31 * C(1,4,K) 

1 132 C(2,4,K)} ¢ 132 

C(4,4,8) (C(4,4,K) - L4l © C(1,4,R) - L42 * €(2,4,K) 
1 - L43 * €(3,4,K)) * Lal 
C(1,5,8) Cil,5,R) * LAI 

C(2,5,K) (C(2,5,R) ~ L22 * C(1,5,K)) * L2I 

C(3,5,K) (C(3,5,8)} L31 Cq1,5,K) 

1 132 © C(2,5,K)) * 231 

C(4,5,8) (C(4,5,8) L461 C(l,5,K} - L42 * C(2,5,K8} 
1 - L43 * €(3,5,K)) * LAI 
U34 * €(4,1,K) 

C24 * C(4,1,R) 

U23 * €(3,1,R) 

Ul4 * €(4,1,K) 

UL3 © €(3,1,K) C(2,2,8) 
U34 * C(4,2,K) 

u2¢ * €(4,2,K) 

U23 * €(3,2,K) 

U1l4 © C(4,2,K) 

V13 * €(3,2,R) C(2,2,K) 
u34 * C(4,3,R) 

uU24 * €(4,3,R) 

023 * €(3,3,K) 

vl4 © C¢4,3,R) 

Ul3 * €(3,3,K) ©(2,3,¥) 
U34 + C(4,4,K) 

Uz4 * C(4,4,F) 

U23 + C(3,4,K) 

ula * C(4,46,K) 

C13 * C(3,4.K} C(2,4,8) 
usa * C(4,5,K) 

24 * C(46,5.K) 

23° €(3,5,F) 


€(3,1,K) © €(3,1,8) 
C(2,1,K) = C(2,1,K) 


1 
C(l,1,&) = CC2,1,K) 
1 


C(3,2,K) c(3,2,KF) 
C(2,2,8) = C(2,2,K) 
1 

C(1,2,K) © Cf1,2,R) 
1 

C(3,3,8) = €(3,3,K) 
C(2,3,8) © €(2,3,R) 
1 

C(1,3,8) C(1,3,K) 
1 

C3,4,8) = €(3,4,R) 
C(2,4,8) = C(2,4,K) 
1 

Ci, 4,R) cq1,4,K) 
1 

C(3,5,K) = C(3,5,K) 
12.5.8) o(2,5,R) 
J 


GO» G8 08 OF OP OD OF OD 
WOnauwawn 


Parte vvr tt rrrpevtrvraree 








—o 


ad al ah nabed 


NAV & Wh: 


aREDEREERRE 








| 
t 


C(l,5,Rp = C(1,5,K) - Ul4 * €(4,5,K) 
1 = U13 * €(3,5,8) - U12 * C(2,5,K) 


3 CONTINUE 
po 6 11 =1,4 
DOO K =2, m1 
9 G6(I1,1,K) = C(I2,2,8) 
DO 6 12 =1,4 
DO6 K -=2, m1 
WH(12,12,1,R) = C(I2,%2+1,K) 
6 CONTINUE 
1000 CONTINUE 
c 
. BACKWARD SUBSTITUTION 
po 71 = IMAX-3,1,-1 
po 7 11 2=1,4 
po 7K = 2, Ml 
GG(I1,1,K) = GG(11,1,K) - BH(I1,1,1,R) 
1 - BE(11,2,1,8) 
2 ~ HA(11,3,1,K) 
3 - HR(I1,4,1,%) 
7 CONTINUE 
RETURN 
END 


SOURCE PROGRAM 


nse2d.f 


GG(1,1+1,KF) 
GG({2,1+1,K) 
GG(3,1+1,K)} 
GG(4,1+1,R) 
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855 : SUBROUTINE MATRX2(IMAX, KMAX) 
656 i COSHH HTERTEEHERAHSSORERRESHAHEKARSSSERAKHSESOHENTARHROTHRERTESSELRIRNAES ES 
857 5 ce * 
“858! ce SUBROUTINE MATRX2 * 
859 1c * 
B60 { CAAeASseea AS eeeRARAHSHERAKRHSESHERKRAASSSHAAEKASST HR ETESESEREPSCOTEXES 
86) : PARAMETER (1X=180,kx=60) 
__ 862 ; COMMON/TRID1/DD(4,4,1%,KX) 
863 , COMMON /TRID2/MM(4,4,IX, KX} 
864, COMMON/TRID3/EE(4,4,1%,KX) 
"865 i COMMON/TRID4/ GG(4,1X,KZ) 
866 | COMMON/SCRAT/A(4,4,1Z) ,HB(4,4,I2,KX) ,C(4,5,IX) 
867 | REAL MM 
BGS | REAL 111,121,131,141,122,132,142,133,143,144 
e659 i 4, LAT, L21, 131,241 
B7 Ihc THIS SUBROUTINE PERFORMS THE BLOCK TRIDIAGONAL MATRIX IVERSION FOR 
5. c AH ENTIRE J=CONSTANT PLANE DURING THE ZETA- SWEEP 
. c 
“874 wa = MAX - 1 
"875 } IM1 = IMAX - 1 
876: po 111 =1,4 
|__ 877 | po1t =2, 10 
878 | ar = 1. / MM(1,1,1,1) 
879: GG(T1,I,1) = GG(I1,I,1) * AI 
8807; WH(I1,1,2,2) = BE(I1,1,1,1) * Al 
881 | BH(I1,2,2,1) = EE(I1,2,1,1) * AI 
82 | WH(I2,3,3,1) = EE(11,3,1,1) * AI 
883 | WH(I1,4,1,1) = BE(11,4,1,1) * AI 
834 | 1 CONTINUE 
c 
886! DO 1000 F = 2 , KMAX - 2 
__887 | Do 5 lls 1,4 
888 | DO 2 Is 2, IM 
889 | C¢11,1,1) = GG(11,1,KR) - DD(12,1,1,K) * GG(1,1,K-1) 
890 | 1 - DD(I1l,2,1,K) * GG(2,1,K-1) 
B91: 2 - DD(11,3,1,K) * GG(3,1,K-1) 
892 | 3 - DD(11,4,1,K) * GG(4,1,K-1) 
"7893: 2 CONTINUE 
894° po5 1281, 4 
"895 | pos I «2, IM . 
"896 A(I1,12,1)= MM(I1,12,1,R) - DD(I1,1,1,K) * BH(1,12,1,K-1) 
$97 | 1 - DD(11,2,1,K) * 8B(2,12,1,R-1) 
! 2 - DD(11,3,1,K) * BH(3,12,1,K-1) 
: 3 - DD(11,4,1,K) * BB(4,12,1,K-1) 
: : CCIL, 1242, 2)= EE(I1,12,1,%) 
901 5 CONTINUE 
902 po3 r= 2, IMl 
903 L1ll = AQ1,21,1) 
904: Lr= 1. * Lil 
_ 905; Ui2 = A(1,2,1) * LIT 
906 . V13 = A(1,3,1) * L1I 
907 | Ul4 = A(1,4,I1) * Lil 
- 908: 121 = A(2,1,T) 
-, 909 | L31 = A(2,1,1) 
910 L141 = A(4,1,1) 
‘91 L22 = A(2.2,1) - L21 * U12 
912: L21 = 1. + 122 
2 G3 | U23 = (A(2,3,1) - L21 * 013) * L2z 
914, U24 © (A(2,4,1) - L21 * Ul4) * L2I 
915, L32 = A(3,2,I) - L31 * U12 
916. L42 = A(4,2,1) - L4l * 012 
917° 133 = A(3,3,1) - L31 * U13 - L32 * 23 
918; L3I= 1. / 133 
919. U34 © (A(3,4,1) - L31 * U14 - L32 * U24) * 131 
920 . 143 = A(4,3,1) - Lal * U13 - L42 * U23 
92) | L44 = A(4,4,1) - L41 * U14 ~ L462 * 0246 - 143 * 034 
922 | Lal = 2. / 144 
923: C(l,1,1) = C(1,1,1) * LLI 
924. C(2,1,1) = (C(2,1,1) - L221 * C¢1,1,3)) * L2r 
. 925 , C(3,1,1) = (€¢(3,1,1) - L31 * €(1,1,1) 
926, 1 - L32 * €(2,1,1)) * L3I 
927 | C(4,1,3) = (C¢4,1,1) - LAL * C(1,1,1) - L42 * €(2,1,T7) 
928 1 - 143 * C(3,1,1)) * Lat 
929 ; C(1,2,1) = €(1,2,1) * Lil 
.. 930 C(2,2,1) = (C(2,2,1) - L21 * C(1,2,1)) * L2I 
"931 €(3,2-1) = (€(3,2,1) - 132 * C(1,2,1) 
932 1 ~ L32 * €C(2,2,1)) * L3I 
933, C(4,2,1) = (C(4,2,1) - L4l * C(1,2,1) ~ L42 © €(2,2,1) 
 934- 1 ~ 143 * €(3,2,1)) * LAT 
935) C(1,3,1) = C(1,3,1) * LIT 
"936: C(2,3,1) = (C(2,3,1) - L2L * €¢1,3,1)) * L2I 
“937 | €(3,3,1) = (C(3,3,3) - 131 * €(1,3,2) 
“938 | 1 - 132 - €(2,3,1)) * L3r 
939 | C(4,3,2) = (C(4,3,1) - L41 * C(1,3,1) - L42 * C(2,3,1) 
940 i - L43 * €(3,3,1)) * Lal 
94! | C(l,4,1) © C¢1,4,1) * Li 
$42, C(2,4,1) = (C(2,4,1) - L212 * C¢1,4,1)) © L27 
~ 943) C(3,4,1T) = (C(3,4,1) - L31 * €(1,4,T) 
9445 1 - L32 * €¢2,4,1)) * 131 
945; CCM, 4,1) = (CC4,4,1) ~ LOL * C(1,4,1) - L42 * €(2,4,1) 
946 | - L43 * €(3,4,1)) * Lat 
947 C(1,3,1) # C(1,5,1) * LIT 
948 C(2,5,1) = (C(2,5,7) - L21 * C¢1,5,1}) * L2r 
- 949 | €(3,5,1) = (C(3,5,1) - L31 * C(1,5,1) 
950 1 - L32 * €(2,5,1)) * L3T 
95! C(4,5,T) = (C(4,5,1) - LAl * C(1,5,1) - L42 * €(2,5,1) 
$2; 1 - L43 * €(3,5,1)) * Lat 
_ 953; C09,1,1) = €(3,1,1) - 034 © €¢4,1,7) 
954. C(2,1,1) = C(2,1,1) - U24 © C(4,1,T) 
_ 955 - 1 - U23 * €¢3,1.1) 
956 Ccl,d,t) = C(1,1,1) - U14 * €(4,1,1) 
957 | 1 ~ U1} * €¢3,1,3) - Ul2 * €(2,1,1) 
958 €¢3,2,1) = €(3,2,1) - U34 © €¢4,2,7) 
959 C(2-2,.1) = €(2,2,1) - 024 * C(4,2,1} 
960. 1 - U23 * €(3,2,7) 
96) €(1,2,3) = €(1,2,1) ~ Ula * €(4,2,1) 
962 — 1 - U13 © €(3,2,1) - 012 * C(2,2.1) 
963 C(3,3,1) = €¢3,3,1) - U34 * €(4,3,7) 
964 €(2,3,1) © €(2,3,1) - U24 * €¢4,3,1) 
- 023 * €(3,3,7) 
- Uld * €¢4,3,7) 
~ U3 © €¢3,3,1) - U12 © C¢2,3,3) 
~ U34 * €(4,4,T) 
- U24 © €(4,4,7) 
- 023 * €(3,4,5) 
- 140° ¢(4,4,1) 
- UL} * €C(3,4,1) - O12 © €(2,4,1) 
- 034 * C(4,5.1) 
- 24 * €(4,5,1) 
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- 023 * C(3,5,3) 

CO1,S5,1) = C(1,5,1) - Ul * C(4,5,1) 

1 = U13 * €(3,5,1) - UL2 * €(2,5,1) 
3 CONTINUE 


Do 6 11 
DO 9! 
GG(11,1,8) 
DO 6 12 
DO 6 I 
BH(T1,12,1,°) 
CONTINUE 
1000 CONTINUE 


BACKWARD SUBSTITUTION 


} 


DO7 EK 
po? 11 
DO 71 
GG(11,1,K) GG(I1,1,K) - BH(I1,1,1,%) * GG(1,1,K+1) 
1 BB(I1,2,1,K}) * GG6(2,1,K+1) 
2 BH(12,3,1,K) * G6(3,1,K+1) 
3 BH(I1,4,1,K) * GG(4,1,K+1) 


aa 
S883: 
Navawn =! A 


i 


i 
sey: 


END 
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SUBROUTINE METRIC 


CecoeenarrcvcasarreeeeaxnaseoeREKese He eRENSHOERERHSSOTERRAO TEEN ENACOES 


ce ; . 
. ce SUBROUTINE METRIC . 
fas] * 


CeeaReccerennececeranneceeuseneeeenRkSeoePRRAEHOURREHSOOEH ERS SOE ERNES 
PARAMETER (1X2180,kx=60) 
COMMON /FIX/OMEGA , EDOT 
COMMON /DGRID/DT, IMAX, KMAX, ITEL, I TEU 
COMMON/GRID1/1( IX, KX) ,2( 1%, KX) 
COMMON /GRID/YACOB( 12, KX) 
COMMON /MTRIX/XIXR(IZ, KE), RIZ(IX, KX), ZETAL( IX, KX), ZETAZ(IZ, KX), 
ARIT(IA, KX), ZETAT(IX, KX} 


: EE 
jsoo5 
aaa OO 
PWM OO RNA 


c 
Cee “YEE SUBROUTINE METRIC COMPUTES THE METRICS IM ALL THE TWO DIRECTIONS AND 
THE GNSTEADY CORFFICIENTS ETAT EFC. 


oils 
a 


DO 1000 K = 1 , EMAX 
DO 1000 I = 1 , IMAX 
ETAU = OMEGA * 2(1,K) 
: YTAU = OMEGA * (-X{I,K))- HDOT 
rt Cees “PRESENT SET UP IS POR FLO® PAST AN AIRFOIL. 


ooICIOo 


IF(1.£Q.1.0R.1.EQ.1IMAX) GO TO 10 


025 | ERI = 65 * (X(141,K)-2(I-1,K)) 
26 BRI = 68 * (2(141,R)-2(1-1,K)) 
“O27 | GO TO 15 
J028 10 IF(I.£Q.IMAX) GO TO 16 
1029 ERI = 1.0 * (X(2,K) — X(1,K)) 
TO3O |} 2X1 = 1.0 * (2(2,K) — 2¢1,K)) 
fo3l: GO TO 15 
‘Jo32 216 EXT = 2.0 © (Z(IMAX,K} - X(IMAX~2,K)) 
033 | ZXI = 1.0 * (2CIMAX,KR} ~ Z2(IMAX~1,K)) 
034, 15 CONTINUE 
O35 ! IF(K.EQ.1.OR.K.EQ.KMAX) GO TO 17 
1036 : RZET © 15 *(K(L,R41)-X(1,K-1)) 
037 i ZZET = .5 *(2(1,R+1)-2(1,K-1)) 
O38 | GO TO 20 
039; 17 IF(K.EQ.RMAX) GO TO 18 
TO40 | AZET © 2. * X(1,2)-1.5 * (1,1) - .5 * ¥(I,3) 
“Yoat | “LUET = 2. * 2(1,2) - 1.5 * 21,1) - -5 * 2(1,3) 


GO TO 20 
18 XZET = 1.5 * RCI, RMAX)-2.* A(T, RMAX-1)4+.5*X(1, KMAX~2) 
ZUET = 1.5 * 2(1,KMAX)-2.* 2(1,KMAX-1)+.5*2(1,KMAX-2) 
20 CONTINUE 
YACOBI © XXI * ZZET ~ XZET * 2X1 
YACOB(I,K) * 1. / YACOBI 
RIX(I,K) = 2ZET » YACOB(I,K) 
XIZ(1,K) © -XZET * YACOB(I,K) 
ATAU = OMEGA * 2(1,K) 
YTAL © - OMEGA * X(1,K)- HDOT 
: RIT(I,K) * - XIACI,R) * ATAU - XLZ(I,K) * YTAU 
$3: 2ETAX(I,K) = ~ZX1I * YACOB(I,K) 
54 ZETAZ(1,K) = XXI * YACOB(I,K) 
5S: ZETAT(1,R) © - ZETAX(I,K) * ATAU - ZETAZ(I,K) * YTAU 
$6: 1000 CONTINUE 
57; RETURN 
$8 END 
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SUBROUTINE DISSIP 

CRAeCRR AA Ceee eKaHe SAKKEASSEERENATHEHAERKESESALERAAGC RT EARTH ESEREAKLE SAS 

ce : . * 

acs SUBROVTINE DISsiP * 
c * 

‘ CeARAAH SCAR AHS SCHERER SHLERTAASEHERERESHOHRENASCEEAEENSCO LEK ESE CEEERED 
PARAMETER (IX=160,kx=60) 
COMMON/FLOW/Q1 (IX, KX) ,Q2(IX, KX) ,Q3(IZ, KX) ,Q4 (1X, KX) 

COMMON /MTRIZ/XIZ (IX, KX), ARIZ (IX, KX), ZETAN(IX, KX), ZETAZ( 1X, KX), 

URTT(IA, RX), ZBETAT(IX, KX} 
COMMON/PERTR/DQL (1X, KX) ,DQ2 (IX, KX) ,DQ3( IX, KX) ,DQ4(1X, KX) 
COMMON/DGRID/DT, IMAX, KMAX, ITEL, ITEU 

COMMON/PAR/GAMMA , REYREF , ALFA, ALFA], REDFRE, AMINF , ALFAI 

COMMON /GRID/Y¥ACOB(IX, KX) . 
COMMON /DAMP/ Wei, WII, WH2X, WHY, WAX, WAY 

COMMON/RADIUS/ SPECX(IX, KX) ,SPECY(1X,KX) 

COMMON/SPEED/ U(IX, KX), V(IX, KX) ,AA(IX, KX) 

DIMENSION P(IX),EPS(1X),DIS1(IX,4},DIS2(IX,4) ,DIS3(IX,4) 
DIMENSION 00(4) 


NOCO OCS! 
NSIS 
FA awn —O 


» (S2ELS SUBROUTINE ADDS TBE FOURTH ORDER DISSIPATION TERMS TO THE 
-MRYQNT BAND SIDE 


IML = IMAX - 1 
IM2 = IMAX - 2 
KMl = KMAX - 1 
KM2 = KMAX - 2 
WOT = WHeDT 
GM1 <= GAMMA - 1. 
DO 20 K = 2, KML 


COMPUTE SWTICHING FUNCITON BASED ON SECOND DERIVAITVE OF PRESSURE 


poldrzI:=1 =, IMAX 
1 P(I) = GM1 * (Q4(1,K) — 0.5°Q1(1,K)*(U(1,K)**2 + V(I,R)**2)) 
po3 r= 2, IM 
Psi = P(I+1l) - 2.*P(1) + P(I-1) 
PS2 = P(I*1]) + 2.*P(I) + P(I-1) 


EPS(I) = ABS(PS1/PS2) 

3 CONTINUE 
EPS(1) = EPS(2) : 
EPS(IMAX) = EPS(IM1) 


SMOOTH OUT PRESSURE COEIFFICIENT 


po4y=2, IM) 

P(I) = AMAX1( EPS(I+1),EPS(1),EPS(I-1) ) 
4 CONTINUE 

P(1) = P(2) 

PCIMAX) = P(IM1) 


Besssossessss 








It 
1108 | 
WW i 
ode 
an, polo rl = 2, IM 
INN | C2 = WH2X * P(T) : 
ED c4 = WW4X - AMIN] (WW4X,C2) 
4a | C22 = C2*WDT * ( SPECX(1,K) + SPECX(I+1,K) ) 
“LS | ca = ~C4*WDT* SPECX(I,R) 
mil wns PY. = O4(I-1,K) + P(I-1) 
milvans =P = Q4(I,K) + P(T) 
Aligbic BPP = Q4CI+1, Ry 4% PEIt1) - 
VIS: DISL¢1,1) = C22" ( Q1lCI+1,R) - Q1(I,K) ) 
420° DIS1(1,2) = €22* ( Q2(1+1,K) - Q2(1,K) } 
marin DIS1(I,3) = €22* ( Q3(I+1,K) - Q3(1,K) ) 
22, DIS1(1,4) = C22* (¢ Q4(1+1,R) - Q4(I,K) ) 
4923 ,¢ DIS1(1,4) = C22" ( BPP - BP ) 
tl24- DIS2¢I,1) = C44 = ¢ Q2CT+2L, K) - 2.°Q2(1,K) + Q1(I-2,K) ) 
12s; DIS2(1,2) = C44 * (¢ Q2(141,K) - 2.*Q2(1,K) + Q2(1-1,K) ) 
1126; DIS2(1,3) = C44 * ( Q3(I4+1,K) ~ 2.*Q3(1,R) + Q3(I-1,K) ) 
27 | DIS2(1,4) = C44 * ( Q4(T41,K) ~ 2.*Q4(1,K) + Q4(I-1,K) } 
128 ye DIS2(1,4) = C44 * ( BPP - 2.°HP + BPM ) 
29: 10 CONTINUE 
Vr30]¢ B.C. TREATMENT 
od FB | QQ(1) = O1(2,K) ~ 01(1,K) 
"V432 | Q0(2} = Q2(2,K) ~ O2(1,K) 
N33 | QQ(3) = O3(2,K) - Q3(1,K) 
A134 | QQ(4) = G4(2,K) ~ O4(1,K) 
135 be QOC4) = 04(2,K) + P(2) - QA(1,R) - PCL) 
1136; DO 15 NA “1,4 
mera C2 = WW2K"P(1) 
1338) C22 = C2 * WOT * { SPECX(1,K) + SPECX(2,K) } 
1139 | DIS2(1,N4) = C22 * QQ(N4) 
“1140 DIS2(1,N4) = 0. 
“V4 15 DIS2(IMAX.N4) = 0. 
“H42}e¢ : 
‘1143 po16ét =1, IM 
i. (h46 DIS3(I,1) = DIS1(I,1) + DIS2(1+1,1) - DIS2{1,1) 
4s DIS3(I,2) = DIS1(1,2) + DIS2(I+1,2) - DIS2(I,2) 
146} DIS3(1,3) = DIS1¢1,3) + DIS2(I4+1,3) - DIS2(1,3) 
ab V42 | DIS3(1,4) = DIS1(1,4) + DIS2(I+1,4) - DIS2(1,4) 
“y148 16 CONTINUE 
‘tHagic 
d 12°- c FILL IN DISSIPATION TERMS 
e 
82 | poi8r = 2, IML 
| $153 | DQ1(1,R) = DIS3({I,2) - DIS3(I-1,1) 
“154 DO2(1,K) = DIS3(I,2) - DIS3(I-1,2) . 
SS DO3(1,K) = DIS3(1,3) - DIS3{I-1,3) 
1136 | DO4(1,K) = DIS3(1,4) - DI$3{1-1,4) 
VIS? | 18 CONTINUE 
VRS | 20 CONTINUE 
Tts9 ic 
“W160 !c K DIRECTION 
Verfc > 
1162 | po 100 r= 2, Imi 
“T1163 ic 
ties ie COMPUTE SWTICHING FUNCTION BASED ON SECOND DERIVATIVE OF PRESSURE 
re 
1166. DO 31 K = 2, KMAX 
| bt67 31 PK} = GM1l * (Q4(1,R) - O.5°Q1(7,K) * (UCI, K)**2 + Viz, KR}**2)) 
1168 | P(1) = ¢ 4.*P(2) - P(3) ) 7 3. 
1169! pO.33 K= 2, RM 
1170, Psi = P(R+1) - 2.°P(K) + PCK-1) 
71 j PS2 © P(K+L) + 2.°P(R) + P(R-1) 
9472) EPS(K) = ABS(PS1/PS2) 
1173! 33 CONTINUE 
Vi74: EPS(1) = EPS(2) 
Ht75: EPS(KMAX) © EPS(KM1) 
WN76i¢ 
fee it SMOOTH OUT PRESSURE COEFFICIENT 
< 
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SOURCE TEXT 


pO 34K = 2, EMI 

P(K) © AMAZ1( EPS(K+1),EPS(K),EPS(K-1) ) 
CONTINUE 

Pcl) = P(2) 

P(KMAX) = P(KM1) 


po 52 Kk 
c2 

ca 

C22 

c44 

—PH 

=P 

PP 
DIS1(K,1) 
DIS1(K,2) 
DIS1(K,3) 
DIS1(K,4) 
DISI(K,4) 
DIS2(K,1) 
DIS2(K,2) 
DIS2(K,3) 
DIS2(K,4) 
DIS2(K,4) 
CONTINUE 
B.C, TREATMENT 

QQ12) ©O1(1,2) - O1(I,1) 

9002) Q2(1,2) - Q2(1,2) 

QQ(3) Q3(1,2) Q3(I,1) 

@ot4) 41,2) ~ O4(1,1) 

QQ14) = O4(1,2) + PYZ) - O4(T,2) ~ PY) 

DO 59. N 1,4 

c2 ww2y © P(1) 

22 C2 * WDT * ( SPECY(I,1) + SPECY(1,2) ) 
DIS1(1,N4) €22 * QQ(N4) 

DIS2(1,N4) 0. 

DIS2{KMAX,N4) = 0. 


2, 

We2y * P(K) 

Way - AMIN1(WN4Y,C2) 

C2 * MDT * ( SPECY(I,K) + SPECY(I,K+1) ) 
-CA*WDT* SPECY(1,K) 

Q4(2,K-1) + P(K-2) 

QO(1,K) + B{K) 

O4(1,K+2} + P(K+1) 

C22 * (¢ QL(I, R41) - Q1(1,K) } 

C22 * ( Q2(1,F+1) - Q2(1-K) ) 

C22 * ( O3¢1,K+1) - Q3(1,K) ) 

€22 Q4(I, R41) - O4(1,K) ) 

€22 wrP ~ #P ) 

ca4 Q1(1,K+1) — 2.°Q1(1,K) + QL(I,K-1) 
cea Q2(1,K+1) - 2.°Q2(1,K) + Q2(1,K-1) 
cea Q3¢1,K+1) - 2.°Q3¢1,K) + Q3(1,K-1) 
C44 Q4(1,R+1) — 2.°Q4(I,K) + Q4(1,K-1) 
C44 SPP — 2.4EP + BPM ) 


— 
ann 


n 


DO 60 K 
DIS3(K,1) 
DIS3(K,2) 
DIS3(K,3) 
DIS3(K,4) 
CONTINUE 


1, RM) 

DIS1(K,1) + DIS2(K+1,1) - D1S2(K,1) 
DIS1(K,2) + DIS2(K+1,2) - DIS2¢K,2} 
DIS1(K,3) + DIS2(K+1,3) - DIS2(K,3) 
DIS1(K,4) + DIS2(R+1,4) - DIS2(K,4) 


anne en 


a 
°o 


FILL IN DISSIPATION TERMS 


NRNNNANNANNANANARNR 
NNN ee eee ew ewe HOO 
VEWN—OBeNAUNEWNKOOM 


aang 


po 65 K 2, KML 
DOQL(I.K) = DQ1(1,K) + DIS3(K,1) 
DQ2(1,K) © DQZ(1,K) + DIS3(K,2) 
0Q3(1,K) = 0Q3(I,K) + DIS3(K,3) 
DO4(I,K) = DOA(I,.R) + DIS3(K,4) 
CONTINUE 

CONTINUE 

RETURN 

END 


DIS3(R-1,1) 
DIS3(K-1,2) 
DIS3(R-1,3) 
DIS3(K-1,4) 











SOURCE PROGRAM DATE 1/20/90 
nse2d.f [TIME 11:46:04 am | 


oe aed 


SUBROUTINE SPECR(ISPEC) 


ae VACRALSTSCARESAOSAERAAESSERATASSHERAKAASSERERSEAOS 
* 
SUBROUTINE SPECR e 
; ce * 
CURA eCCeANA Te ec ERKKSSSOURENASTHTUTENTOTETHESSCO ERT ERGO OTERERACOERER 
PARAMETER (12=180,kx=60) 
COMMON /PAR/GAMMA , REYREF ALFA, ALFA], REDFRE , AMINF , ALFAI 
COMMON/DGR:D/DT, IMAX, KMAX, I TEL, ITEU 
COMMON/MTRIX/ XIX( IX, KX), AIZ( IR, KX), ZETAR(IX, KX), ZETAZ (1X, KX) 
1, 2IT( IX, RX), BETAT (IX, KX) 
COMMON /FLOK/Q1 (IX, KX) ,Q2(12, KX), Q3( 1A. KX) -Q4 (1X, KX) 
COMMON/RADIUS/ SPECX(1X,KX) ,SPECY (1X, Kx) 
COMMON/SPEED/ U(IX, KX), V(IX, KX} AAC IX, EX) 
COMMON/GRID/YACOB(1X, KX) 


“QETS SUBROUTINE COMPUTE THE SPECTRAL RADIUS FOR SCALING TRE 
RXPLICIT AND IMPLICIT DISSIPATIONS 


GAMMA * { GAMMA - 1. ) 

= 1, KMAX 

= 1, IMAX 

Q2(1,K) / QLiI,K) 

V(TOR) © Q3(1.K) / QU(I,K) 

AA(I, K)= GAM * ( Q4(1,K)/Q1(1,K) ~ 0.5*(U(I,K)**2 + V(I,K)**2) ) 
IF(AA(I,K).LT.0.) PRINT*, "NEGATIVE A*A = ',AA(T,K),' AT ', TK 
CONTINUE 


K 1 
1 1 
* ‘ 


Sey 
Ne: 


pete 
Poe 


in 
me 


COMPUTE IMPLICIT DISSIPATION SCALING 
SPECK » SPECTRAL RADIUS FOR XI-DIRECTION 
SPECY = SPECTRAL RADIUS FOR ZETA-DIRECTION 


IF(ISPEC.EQ.1) THEN 
pO 20 K = 1 , KMAX 
po 20 1 1 =, MAX 
SPECK(I,R) = 1. / YACOB(I,R} 
SPECY(I,K} = 1. / YACOB(I,K) 
CONTINUE 
ELSEIF(ISPEC.EQ.2) THEN 
po 30 kK = 1 , 
po 30 1 * , MAX 
UCON = UCT, K)*XIX(I,K) + VI, R)*xIZ(1, KR) + XIT(I,K) 
CON «= U(1,R)*ZETAX(1,K) + V(I,R)*ZETAZ(1,K) + ZETAT(I,K) 
SPECX(I,%) = ABS(UCON) / YACOB(I,K) 
SPECY(1,K) = ABS(VCON) / YACOB(I,K) 
CONTINUE 
ELSEIF(ISPEC.E 
po 40 K = 1 , KMAX 
po 40 r = 1 , IMAX 
UCON = U(I,K)*XIX(1,R) * VC, R)*XIZ(1,K) + XIT(1,R) 
= 
= 
a 
. 


Oonauwawny—: 


Seeey 
tr} epat 2} 


Q.3) THEN 


VCON UCI, K)*ZETAX(I,K) + V(1,K)*ZETAZ(1,K) + ZETAT(I,K) 
X12 MIX(I,K)**2 + REZ, R)e*2 
ZETA2 ZETAX(1,K)**2 + ZETAZ(1,K)**2 
SPECX(I,F) (ABS(UCON) + SQRT(AA(I,K)*XI2) ) 7 YACOB(I,R) 
SPECY(1,K) {ABS(VCON) + SQRT(AA(I,K)*2ETA2} ) / YACOB(1,K) 
CONTINUE 

ENDIF 

RETURN 

END 





SOURCE PROGRAM DATE 1/20/90 
nsezd.f [ie 11:46:06 am | 


SUBROUTINE EIPBC(CL) 
Crceseusnree eteuae PRA TAKS SCE RRASSCHCERAHSSHSHARAHESEHERKHRASCERNKEOOE 
c* t 
jee SUSROUTINE EXPBC * 
ce t 
CRUSH SSKER AAAS SCEERRAKEHAEARTASHHETRAASCERAKAARSSSARERASSHUERERESCOERAE 
PARAMETER (IX=160,kx*60) 
COMMCN/SURF/PSUR (1X) 
COMMON/GRIDL/X (IX, KX), 2(1X, KX) 
COMMON/PAR/GAMMA, REYREF , ALFA, ALFA] ,REDFRE , AMINF 
COMMON/DGRID/DT , IMAX , KMAX , ITEL, ITED 
COMMON /GRID/YACOB (1X, EX) 
COMMON/DAMP/WH, WWI, WW2X , WW2Y , WAX, WAY 
COMMON /MTRIX/XIX (IX, KX) ,XIZ(IX, KX), ZETAR(IX, KX), ZETAZ(IX, KX), 
1RIT(IR, KX), ZETAT( IX, KX) 
COMMON /FLOW/Q1 (1X, RX) ,Q2(IX, KX) ,Q3( IK, KX) ,Q4(IX, KZ} 
COMMON/SPEED/ U(1X,KX),V(1X, KX) ,AA(IX, KX) 
COMMON/INITI/UINF , VINF 
COMMON /BC L0G /CIRCOR 
COMMON/LOGIC/RSTRT , PITCE , RAMP 
LOGICAL RSTRT, PITCE, RAMP 
LOGICAL CIRCOR 
DIMENSION €1(4) 
DIMENSION A(2,2),RHS(2) 
DATA P1/3.1415927/ 
DATA CHORL/1./ 


‘OO Sazy y 


BNAVAUNR 


Soe vous wn, 


ome 


| 
gee 


~ ai 


i 


a 


GAMI = 1. / GAMMA 
GAM™ == GAMMA ~ 1. 
GAMMI = 1. / GAMM 


INVISCID AND VISCOUS B.C. ON SOLID WALL 


aaa 


DO 1I= ITEL , ITEU 

Ke« 3 

C1l(1) = XITCI,R) 

C1(2) = XIX(I,K) 

C1¢3) = ¥IZ(1,R) 

UCON3 (Q2¢1,K)*C1(2)+Q3(1,K)*C1(3)) 
1/01(1,K) 

K= 2 

Cl¢1y = XIT(ILR) 

C1(2) = XIX(I.K) 

C1(3) = XIZ(1,K) 

UCON2 = (Q2(1,K)*C1(2)4+Q3(1,K)*C1(3)) 
1/01(1,K) 

RHS(1) * 2. * UCON2 ~ UCON3 - XIT(I,1) 
POR VISCOLS FLOKS SET UCON TO ZERO ALSO 
IF(REYREF GT.0.) RHS{1) = - XIT(I,21) 
A(1,1) * XIX(1,1) 

A(1,2) = X12(1,1) 

A(2,1) = 2ETAX(I,1) 

A(2,2) = 2ETAZ(12,1) 

RHS(2) = - ZETAT(1,1) 

TEMP) A(1,1) 

TEMP2 = A(1,2) 

TEMP3 = A(2,1) 

TEMP4 = A(2,2) 

DEN = 1. /(TEMP1 * TEMP4 - TEMP2 * TEMP3) 
A(L,1) * A(2,2) * DEN 

A(1,2) = - TEMP2 * DEN 


eee nanan SOE aie ales Sneaanee Slr mae? 


N= COONAVEWN —OOSNAUAWN= 


a ae a pe ee fe re 
baDS SPWWUW WWW WWW NRA IN NIRS 
w 


if 
‘ 
‘ 
: 


A(2,1) = - TEMP3 * DEN 

A(2,2) = TEMP] * DEN 

Q1(I,1) = 2. * QI(1,2) - O1¢1,3) 

Q2(1,1) = QLC(I,1)*¢A¢L,1)*RHS(1)*A(1,2)*RHS(2)) 
Q3(1,1) = QU(I,1) *(A(2,1)*RBS(1)+A(2,2)*RHS(2)) 
CONTINUE 

DO 10 IeITEL ,ITEV 

P2=GAMM* (04(1,2)-0.5°Q1 (1,2) *(U(I,2)**24¥(1,2)**2)) 
P3=GAMM*(04(1,3)-0.5*Q1(1,3)*(U(1,3) **24V(1,3)**2)) 
Pl=(4.*P2-P3)/3. 

PSUR(I)= (GAMMA*P1-1.)/(.7*AMINF**2) 
Q4(1,1)"P1/GAMM+0.5*Q1(1,1)*(U(1,1)**24+¥(1,1)**2) 


5 
5 
5 
5 
5 
5 
5 
5 
5 
5 


350 
35) 
352 
353 
1354 
355 | 
356 
357; 
358 
359 
360. 
361: 


ae 
a 


WOnaAvVawn—OOa 


FAR FIELD BOUNDARY CONDITION ONLY FOR STEADY FLOW 


circ = 0 
IF(PITCR.CR.RAMP) GO TO 999 
IF(AMINF.GT.1.) GO To 65 


a ce ee ce ee ee ee ee oe 


Dds bdr bee de bed, 4 dp aded abadebebabeb6 


CIRCULATION CORRECTION AT TBE FAR FIELD IS BASED ON POTENTIAL 
VORTEX 


PAE AE A AP OT PA AT ad Sad 


WOCCOHO OOO SPSSSPR wae 


~NNNNNNNNNO 


BETA = SORT( 1. - AMINF**2 ) 
IF(CIRCOR) CIRC = 0.25 * CHORD * CL * BETA * AMINF / PI 


COSAL = COS(ALFA) 

SINAL = SIN(ALFA) 

AINF #1. 

BINF = GAMMI + 0.5 * AMINF**2 

CIRCQ = CIRCULATION CORRECTION QUANTITY 


ne ne ie eee ee ee ne ee ee we ee ee 


aphbbR RS 


SSSSR rau 


K = EMA 

Do 60 I = 2 , IMAX-1 

XLOC = X(1,KR) - IREF 

ZLOC = %(1,K) 

RADIUS = SORT( XLOC**2 + ZLOC*#2 } 

ANGLE = ATAN2(ZLOC,2L0C) 

CIRCQ * CIRC / ( RADIUS * ( 1. - (AMINF*SIN(ANGLE-ALFA)}**2 ) } 
UF = UINF + CIRCQ * SIN(ANGLE) 

VF = VINF - CIRCQ * COS(ANGLE) 

AFSQ = GAMM * ( BINF - 0.5*( UF*#*2 + VFe#2 ) ) 
ar = SORT(AFSQ) 


NOWREFLECTING B.C. BASE ON 1-D RIEMANN INVARIANTS 
BETXN, SETZN = NORMALIZED NORMAL VECTOR 


ANOR = 1. / SQRT( ZETAX(I,R)°*2 + ZETAZ(I,K}**2 ) 
ZETXN = ZETAX(I,R) * ANOR 
ZETZN = ZETAZ(1,R) * ANOP 


SSSSSSSSES 
DAWN HK$ OCOBNAW WN KOU ONGHUARWN 


CHECK FOR INFLOW OR OUTFLOW 

FOR INFIOK- Rl, VELT, ENTROY ARE SPECIFIED AS FREE STREAM VALUES 
R2 1S EXTRAPOLATED 

FOR OUTFLO®: R1 15 SPECIFIED AS FREE STREAM VALUE 
R2, VELT, ENTROPY ARE EXTRAPOLATED 


AANANnAN 


RBOEXT = Q1(3,R-1) 
RAO? =) Ot(1.%-1) 










SOURCE PROCRAM 












17 VEXT «= 7 2(1,R-1) * REOI 
are | VEXT + ¥3(1,R-2} * RHOI 
“1419 | EEXT = 94(1,K-i) 
#0, PEXT = GAMM * ( EEXT - 0.S*RHOEXT*®( UEXT**2 + VEXT**2 ) } 
;c 
Waazic SET RIEMANN INVARIANTS &1, AND R2 
a is VELN = NORMAL VELOCITY, VELT = TANGENTIAL VELOCITY 
fe 
4355 Rl = ZETXN * UF + ZETZN * VF - 2. * AF * GAMMI 
| £426 | R2 = 2ETUN * VEXT + ZETZN * VERT + 2. * SQRT( GAMMA © PEXT / 
183 1 RHOEXT } * GAMMI 
ic 
4429 | VELN = ( R1 + R2) * 0.5 
1430 SPSQ = ( R2- Rl) * GAMM * 0.25 
lay A2 = SPSQ**2 
1432 [Cc 
1433 c SET OTHER FIXED OR EXTRAPOLATED VARIABLES 
c 
1435 | IF(VELN.LE.0.) THEN 
1436 | VELT © ZETZN * UF ~- ZETXN * VF 
4 432 | ENTRO = GAMMA 
4 ! ELSE 
aie VELT = ZETZN * UEXT - 2ETAN * VEXT 
bord | ENTRO = RHOEXT**GAMMA / PEXT 
i ENDIF 
442 is 
443 bc NOW COMPUTE FLOW VARIABLES 
f laaaic 
V4a4s |; U(I,K) = ZETXN * VELN + ZETZN * VELT 
1446; V(I,K) © 2ZETZN * VFLN - ZETXN * VELT 
1447 | Q1(I,K) = ( A2 * ENTRO * GAMI ) ** GAMMI 
1448: PRESS = A2 * Q1(1,K) * GAMI 
1449 | Q2(1,K) = Q1(1,K) * U¢I,R) 
Y $450: Q3¢1,K) = QU(I,K) * V(I,K) 
1450; Q4(1,K) = PRESS*GAMMI + 0.5 * Q1(I,K} * ( UCI,K)* ¢ + V(1,K)**2 ) 
1452: 60 CONTINUE 
1453: GO TO 67 
eee 65 CONTINUE 
eo 
1456 6c B.C. FOR SUPERSONIC FLOW 
1488S 
4 K = KMAX 
1459. DO 66 I = 2 =, IMAX-1 
1460 RHO] © 1. / Q1CI.R) 
1461. UCI,R) = Q2(1,K) * REOI 
1462 . VCILR) = Q3(1,K) * RAOI 
, 1463 | ANOR = 1. / SORT( ZETAX(1,K)**2 + ZETAZ(I,&)**2 ) 
1464 | ZETXN = ZETAX(I,K) * ANOR 
1465 ZETZN = ZETAZ(1,K) * ANOR 
1466 | VELN-* ZETXN * U(I,K) + ZETZN * V(I,R) 
1467 : TE(VELN.GE.0.) THEN 
1468 | QL(I,K) = Q1(1,-1) 
1469 Q2(1,K) = Q2(1,R-1) 
1470 Q3(1,K) = Q3¢1,K-1} 
1470 ; O4(I,K) = Q4(I, K-12) 
1472. ENDIF 
pa 66 CONTINUE 
:< 
1475 | 67 CONTINUE 
1476 °c 
te ae OUTFLOW B.C. AT THE DOWNSTREAM OF C-GRID ONLY FOR INVISCID FLOW 
| Cc 
1479 | IF(REYREF.LT.0.) THEN 
1480 | I= 
4481 - IPl = 1 
ered SIGN * -]. 
ic 
1484 ¢ CHECK FOR SUPERSONIC FLOW 
1485 ;¢ 
1486 | IF(AMINE.GT.1.) GO TO 75 
1487 6c 
1488 | 72 CONTINUE 
ae DO 74 K = 1 = , KMAX 
te 
ae {c CORRECT FREE STREAM VALUES WITH CIRCULATION CORRECTION 
c 
4 
1493 | XLOC = 3(1,K) - XREF 
1494 zoe = 2¢1,K) 
1495 | RADIUS = SQORT( XLOC**2 + ZLOC#*2 } 
1496 | ANGLE = ATAN2(ZLOC,XLOC) 
1497 , CIRCQ = CIRC / ( RADIUS * ( 1. - (AMINF*SIN(ANGLE-ALFA))**2 ) ) 
$498 | UF © UINF + CIRCQ * SIN(ANGLE) 
1499 | vF = VINF - CIRCQ * COS(ANGLE) 
1500 | AFSQ = GAMM * ( BINF - 0.5°( UF**2 + VF**2 ) } 
: eat | AF = SORT(AFSQ) 
ie 
4 203 c RINE, RIZH = NORMALIZED HORIZONTAL VECTOR 
¢ 
($65 ANOR = 1. / SORT( RIX(I,K)**2 + XIZ(I,R)**2 ) 
1506 | RIRB == RIX(I,K) * ANOR 
eee RIZB o* XIZ(2,K) * ANOR 
ic 
1$09 RAOEXT = Q1(1+IP1,K) 
1510; RROl = 1. / QL(1+IP1,R) 
iio UVEXT == Q2(141P1,K) 
1$t2! VEXT © Q3(1+IP1,K) 
$513, EEXT == Q4(1+IP1-K) 
rie PEXT © GAMM * ( EEXT ~ 0.5*REOEXT*( UEXT**2 + VEXT**2 ) } 
ic 
1318 ic BET RIEMANN INVARIANTS R1, AND R2 
jc 
1518, Rl = XIMH * UF + XIZH * VF - SIGN * 2. * AF * GAMMI 
1519, R2 =< AIXH * UEXT + XIZB * VEXT + SIGN * 2. * SORT( GAMMA * PEXT / 
lear | 1 RBOEXT ) * GANMI 
¢ 
1§22_ VELN = ¢ R1 * R2) 2 0.5 
1823, SPSQ = SIGN * ( R2 - Ri ) * GAMM * 0.25 
bese A2 = SPSQ°*2 
ic 
ey ic SET OTBER FILED OR EXTRAPOLATED VARIABLES 
Cc 
128° IF(SIGN*VELN.LE.0.) THEN 
1829 VELT = -XIZR * UF + XIXH * VF 
1530 — ENTRO = GAMMA 
1531. ELSE 
1532 | VELT © -XIZR * UEXT + XIXA * VEXT 
($33 - ENTRO = RBOEXT**GAMMA / PEXT 
ie ENDIF 
1236 COMPUTE FLOK VARIABLES 





SOURCE PROGRAM DATE 1/20/90 ee _zeren Mg | 


nse2d.f TIME | TIME 11:46:04 am | 


eae a a a ae 
SOURCE TEXT 


U(I,K) = XIXB * VELN - XI2H * VELT 
V(1,K) = X12B * VELN + XIXH * VELT 
Qk(I,K) = ( A2 * ENTRO * GAMI ) ** GAMMI 
PRESS = A2 * Q1(1,K)} * GAMI 

Q2(1,K) = Q1LCI,K) * U(I,K) 

Q3(I,K) = QLU(I,K) * V(1,K) 

Q4(1,R) = PRESS*GAMMI + 0.5 * QI(I,K) * ¢ UCI,R)**2 + V(I,R)**2 ) 
CONTINUE 

IF(I.EQ.IMAX) GO TO 79 

I = IMAX 

Ipl = -1 

SIGN = 1. 

GO TO 72 

CONTINUE 

GO TO 89 


tod 


Parvad ak ad ad ah vad adalah aitad ah ai ak at alas 


SASSAae 


B.C. POR SUPERSONIC FLOW 


aes Bs re ca, 


CONTINUE 
pO 84 KF = 1 = , KMAX 


VVVAIvAW a 
WODNAMNAWN =O” 


RBCI «1. / O1(I,K) 

UI, KR) = Q2¢1,K) * RHO! 

V(1,R) = Q3(1,R) * RHOI 

ANOR 1. / SQRT( RIX(1,K)**2 + XIZ(I,R}**2 ) 
X18 XIX(I,K) © ANOR 

XIZB XIZ(I,R) * ANOR 

VELN SIGN © ( ZIXB * (I,K) + XI2B * V(I,R) ) 


Vw 
oneona 


a 


IF(VEIN.GE.0.) THEN 
“Q1(1,K) = QL(I+IP1,K) 
Q2(I,K) = Q2(I+IP1,K) 
Q3(1.K) = O3(1+IP1,K) 
Q4(1,K) = Q4(1+IP1,K) 

ENDIF 

CONTINUE 

IF(I.EQ.1IMAX) GO TO 89 

I =» IMAX 

IPl = ~1. 

SIGN = 1. 

GO TO 75 

CONTINUE 


DOUNSTREAM B.C. FOR VISCOUS FLOW 
ELSEIF(REYREF.GT.0.) TREN 
OUTFLOW B.C. FOR VISCOUS FLOW 


RHO, U AND V ARE SET TO THE VALUES OF TRE NEXT INTERIOR POINTS 
PRESSURE 1S EXTRAPOLATED AND THEN COMPUTE ENERGY (Q4) 


NANANA AAN 


pO 100 K = 1 , KMAX 
Ie l 

QUC TK) = QL(I+1,K) 
Q2(1,K) = Q2(1+1,K) 
QUI, K) = Q3(I+1,8) 


PEXT = GAMM * ( Q4(T+1,R) - 0.5 * (Q2(14*1,K)"*2 + Q3(I+1,K)**2) 
1 y QU(1*1,K) ) 

O4(1,F) = PEXT“GAMM + © 5 * ¢ Q2(1,K)**2 + Q3(1,K)**2 ) / QL(I,K) 
I = IMAX 
Ql(I.K) = 

Q2¢1-K) = Q2(1-1.K) 


Ql¢I-1.K) 


Q3(1,R) Q3(1-1,K) 

PEXT GAMM © ( Q4(I-2,R) - 0.5 © (Q2(T-1,K)%*2 + Q3(I-1,K)°*2) 
1 f# QL(I-1,K) ) 

Q4CT.R) = PEXT/GAMM + 0.5 * ( Q2(1,K)**2 + Q3(1,K)**2 ) / QLCI,R) 
CONTINUE 


ENDIF 
AVERAGE FLON VARIABLES ACROSS WAKE CUT (FOR C-GRID) 


g 


WaWN—-OWOen IVAaAWN—O”O 


CONTINUE 
bo 120 7 1, #%ITEL- 1 
Iv = IMAX + - JI 
QIAVG = 0.5 { Ql(I,2) + QU¢IU,2) } 
Q2ave = 0.5 (€ Q2¢1,2) + Q2qtl,2) ) 
Q3AVG = O.5 * ¢ QU(1,2) + Q3(IU,2) ) 
QlcT.1) = QLAVG 
Q1(1U,1) © QLAVG 
Q2(3,1) = Q2ave 
Q2(1U,1) = Q2AVG 
Q3(1,2) = Q3AVG 
Q3(TU,1) = Q3IAVG 
PSLOW = GAMM * ( Q4(1,2} - 0.5 * ( Q2¢1,2}**2 + O3(1,2)**2 5 / 
1 Q1¢1.2) ) 
PSUPP = GAMM * (¢ QO4(IU,2) - 0.5 © ¢€ Q2¢1N.2)-*24Q3(1U,2)**2 ) / 
1 Q1(IU,2) ) 
PSAVG = 0 5 * ( PSLOW + PSUPP } 
Q4AVG © PSAVG » GAMM + 0.5 * ( Q2(1,1)**2 + Q3(1,1)**2 ) / 
1 Q1(1,1) 
Q4(I.1) = Q4avG 
Q4(1U.1) = Q4AVG 
120 CONTINUE 
RETURN 
END 
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SOURCE PROCRAM DATE 1/20/90 
nsezd.f : 
= 5 : 






1634 SUBROUTINE STRESS(ITN) 

1635 | CrcccnunasecesacancecranmirevterRtecerTEHAeeeERRNSDOeRreeeesReKerer 
1636 ce * 
1637 { ce SUBROUTINE STRESS . 
1638 ce * 
1639 i Caen aeecceannsevecernsaeetr AACE RET RHCOCERAERECCERRROOCEEERRHO CETERA 
640: PARAMETER (IX*1°0,kx*60) 

COMMON /FLOW/O1 (1X, KX) ,Q2(IX, KX) ,Q3(1Z, KX) ,Q4(.2, KX) 
COMMON/DGRID/DT, IMAX, KMAX, ITEL, ITEU 
COMMON/GR.D1/X (12, KX), 2(12, FX) 

COMMON /PARK /GAMMA, REYREF ALFA, ALFAL,REDFRE , AMINE ,ALFAI 

: COMMON/SPEED/ U(IZ, KX) ,V(2%, 8X) -AA(14, RX) 

1646 ; COMMON/PERTR/DOQL (1X, KX), DQ2 (IX, EX), DQ3( 1%, KX) ,DQ4(.%, KX) 


1647 ; COMMON/MUTUR/CMU (IX, RX) 

1643, DIMENSION RB1(1X)},RH2(1X),RB3( 12) ,RBA( 1X) 

1649 | COMMON/LOGIC/RSTRT , PITCH, RAMP 

1650 | LOGICAL RSTRT, PITCH, RAMP 

is }- c THIS SUBROUTINE ADDS VISCOUS TERMS TO THE RIGHT BAND SIDE 

c 

Vas c CMUL = LAMINAR OR MOLECULAR VISCOSITY 
18 c 

1655 ; CMUL = AMI. > / REYREF 

1656! ¢ CALL EDDY (CHUL) 

1657: RWAXML © KMAX - 1 

Voces IMAZM1 + IMAX - 1 

PR 1. 

1660 | c 

ies te COMPUTE TXX, TRY AND VISCOUS DISSIPATION AT I~ 1 / 2 

i] rc 

1663. DO 30 K = 2 , RMAXM) 

1664 | BDO 20 1 = 2, IMAX 

1665 | UXT = UCI, K) - UyI-1,K) 

1666 ; VXI © V(ILK} - V(I-1,K) 

1667 : AXI = AA(i,R) - AA(I-1,R) 

1668 | U2ET= .25: (UCT, R+L)-UC1,F veI-1,R+1)-U¢1I-1,R-1)) 

1669 | VZET= 125° (VCE, R+L)-V(1,R-2) 0. 2-1, 841 )-V(I-1,K-1)) 
1670. AZET= .25°(AA(1,K*1)~AA(1,K-1)+AA(I-1,K+1)-AA(I-1,K-1))} 
1671. ANI = X(T. K) - N(I-1,K) 

1672, ZX1 = Z1,K) - 2¢1-1,8) 

1673 RZET= .25 © (X(T, K+h)-X(1,R-1)+2(1-1,K+1)-X(1-1,R-1)) 
1674 VLET= 28 + (21, Kl )-2¢2,K-1)4+2(1-1,Be2h)-Z(1-1,K-1)) 
4675. YAC = AXI1 * Z2ET - 2X1 * XCET 

1676, Yac = 1. yac 

1677. XIX = 2ZET * YAC 

1678 | ZETAX= - 2XI * YAC 
| 1679. RIZ = -XTET * YAC ‘ 
1680 | ZETAZ= AI * YAC 

1681 CNM = ( 5 * (CMU(I,K) * CMU(I-1,K)) + CMUL ) * DT 

1682 UX = UXI * XIX + UZET * ZETAX 

1683 Vx = VE) * XIX + VZET * ZETAX 

1684. AX = AX] © XIX + AZET * ZETAX 

1685 UZ = UXT * X12 + UCZET © ZETAZ 

1686 | V2 o= VAI * X12 + ET * ZETAZ 

5687. AZ = ARI] * X12 + AZET * 2ETAZ i 
1688 | TAX = -(-4. * UX + 2. * VZ) * CNM / 3. 

1689 TRY © CNM © (U2 + VX) 

1690 | TYY = -ChY / 3. © (-4. 8 V2 + 2. © UX) 

1693 : RE © CCUT,F)+UC I-31, K) pS *TAHe Ev (1, K)+V{1-1,8))*TXY) 0.5 
1692 1 * CNY / PR/(GAMMA ~ 1.) © AX 

1693 SO = (CLT RYSUC T=, RY) STXY4 (VC I, K)+V(I-1,K)) *TYY)*0.5 


+ CNY / PR / (GAMMA - 1.3 © AZ 


$ 

a 
anaananan 

~ 


1695 DEBUG 
1696 | TURN OFF ENRGY DISSIPATION AND DIFFUSION 
1697 | R4 = 0. 
, 1698 - S4 = 0. 
1699 ; RBL(1) = 0. 
1700 ; RA2(1) © (XIX * TXX + NIZ * TRY) / YAC 
1701; RE3(1) = (XIX * TRY * XIZ * TYY) / YAC 
1702 20 RH4(1) = (XIX * R4 + XIZ * $4) / YAC 
1703 . DO 30 1 = 2, IMAXM 
1704 | DQI(T,R) = DQL(I,K) + RB1(1+1) - RB1(T) 
1705 ; DQ2(1, KR) = DQ2(1,K) + RH2(1*1) - RB2(1) 
1706 DQ3(1,8) = DQ3(I.K) + RB3(I+1) - RB3(1) 
1707 | DO4(T, KX} = DOG(1,K) + RB@{I-1) - RB4(]) 
1708 | 30 CONTINUE 
1709 ;¢ IN THE 2 DIRECTION 
1710: po 70 I = 2, IMAXM1 
Wee: DO 60 K = 2, KMAX 
712. URI = 625 © (UCT4+1, R)-UCI-1,K)+U(T+1,R-2)~-U(1-1,8-1)) 
713, VXI = 6.250% (VCTSL,RY-VQI-1 F}eVCTe1, R-1)-V(1-1,8-1)) 
74, AXT @ 125 * (AA(I#1,R)-AA(I-1,K)+AA(141,R-1)-AA(I-1,R-1)) 
1705 | EMT = 125 * (M(142,K)-X(I-1,R)+3(141,8-1)-X(1-1,F-1)) 
| 1716) ZAI = 625% (2CT+) &)Y-2C1-2,R)+2(141,K-1}-2¢1-1,8-1)) 
Viz UZET = U¢I,R) - (I, R-1) 
718) VZET = V(I,K) - V¢I,R-1) 
mele ALET = AA(I,R) - AACT,R-2) 
1720 | RLET = X(T, KR) - ¥(3,K-1) 
Wet. 2LET = 210,K) - UCI, Kk-1) 
1722 | YAC © XX1 * ZZET - 2% * RZET 
1723. yaAC = 1. / YAC 
1724) RIX + LZZET * YAC 
1725) ZETAX= - ZX1 * YAC 
1726! RIZ © -X2FT * YAC 
1727; ZETAZ© XXt * YAC 
1728; CMM = ¢ 25 * (CMU(I,R) * CMU(T,R-1)) + CMUL) * DT 
1729, UR = UR1 * RIX + UZET * ZETAX 
_ 1730 | VE © VXI * XIX + VZET * ZETAX 
W731: AX © AXI ¢ XIX + AZET © ZETAX 
1732 | U2 = UXI * X12 + VIET * 7ETAZ 
1733 | 4% © VAI * XIZ 4 VZET © ZETAZ 
1734; AZ © ARI * XIZ + AZET * ZETAZ 
$735 | TAX = -(-4. 9° UK 4 2. 6 V2) ¢ CNM 3. ‘ 
(736 | TAY = CNY * (UZ + VX) 
3737: TYY © -CNM / 3. * (~@. © VZ + 2. © UX) 
1738: ¢ RE | CKUCT RpsUCT R-L))* TRV CTL ROC, R-1) °TXY) 20.5 
1739 ¢ 1 + CNM / PR/(GAMMA - 1.) * AK 
1740 (6c SE = (CUT RY*UCT , E-2 5) *TRY*(V(1, K)4V(1,R-1} )*TYY)*0.5 
W4lie 1 + CAM / PR / (GAMMA - 1.) * AZ 
1742 | Ra = 0 
$743 | sq=0 
1744) RAL(K) = 0. 
1745 RA2(K) = (ZETAX * TRX + ZETAZ * TRY) / YAC 
1746 RND(K) «= (ZETAX * TRY + ZETAZ * TYY) /% YAC 
1747 | 6O RB4(K) = (ZETAX * R4 + ZETAZ * S4) / YAC 
$748) DO 70 KR = 2 , RMAX™1 
9749, DQ1(1, KR) = DQICI,K) + RBI(K+1) - RB1(R) 
1750 DO2(T, Ry) = DQI(1,K) + RH2CK+i) - RH2(F) 
750 DO':1,R) = DQI(I,R) * RB3(K-1) - RB3(K} 
1782 DQ4(i. KR) = DQOG(I,K) + RB4(K*1) - RB4A(K) 


1753" 70 CONTINUE 





{ SOURCE PROGRAM DATE 1/20/90 
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SOURCE PROCRAM 


nse2d.f 


DATE _-1/20/90 | PACES 









1757 | SUBROUTINE LOAD(CL,CDP,CDF,CM ALFAS , XREF) 
1758 | cocecceunccceunsasevereraucecenanseecensseeeeersnaeeutsEEeeyeenhaseee 
1759 | ce * 
“1760 | c+ SUBROUTINE 2OAD e 
T76) ie . 


V7G2 | CeenravceceareacccanaseceeneneseeeaRRELOTERAROSOCeRTENSCTE ENESCO TEES 





“1763 ' PARAMETER (1X=180,kx=60) 
T7664 ; COMMON /GRIDIL/X(IX, KX), ¥ (1X, KX) 
76S : COMMON / SKINCF/CF (1%) 
1766 : COMMON /DGRID/DT, IMAX, KMAX, ITEL , ITEU 
7 367 | COMMON /SURF/PSUR (IX) 
ic 


W769 ic THIS SUBROUTINE COMPUTES YEE INVISCLD CONTRIBUTIONS 
Wroic TO LOADS ON THE ASRFOIL SURFACE 


“7c 

17721 CL = 0. 

“1773 | cD = 0. 

V774) CDF = 0. 

"W775 cM = 0. 

1776 } DO 300 I = ITEL , ITEU ~ 1 

1277: DX = 2141.1) - X(1,1) 

“1778; «300 CDF = CDF + ( CE(I) + CF(I+1) ) * 0.5 * Dx 
"1779! DO 400 Y = ITEL , ITEU ~ 1 

1780 | AL = 65 * (N(1,2)4X(141,1)) 

T78r ; YL = £8 * (¥(1,1)+¥(1+1,1)) 

1782 : DX = X(I+1,1) - ¥(1,1) 

1783: DY = ¥(I+1,1) - Y(I,1) 

1784) CPA = PSUR(1I+1) * .5 + PSUR(I) * .5 
_1785 | DCL = CPA * (-DX) 

“1786 | pcp = CPA * DY 

“1787 | CL = CL + DCL 

1788 | ch = CD + DCD 

1789! 400 CM = CM + DCD * YL - DCL * ( ZL - XREF ) 
"1790 ic 

W791 | DCL = CL * COS(ALFAS) - CD * SIN(ALFAS) 
1792: CDP = CL * SIN(ALFAS) + CD * COS(ALFAS} 
1793 | CL = DCL 

“1794: RETURN 


1795 END 





| 





SOURCE PROCRAM 
nse2d.f 
2. Sa z 


SUBROUTINE WRAP(II,J3,XSING,YSING ,XP,YP,SO0,A0,YO) 

CheeeeaRKeaseha wR Kee eeAAHSHORARAHSSSUTAAHSSUTEHTRSSSSERARHTEETHRAESEE 

* 

SUBROUTINE WRAP ‘ 

* 

COCR SHSEET RES SEAEREHHOEERAROOSHSE REAP HELEAESHSCEETEASHERERETEOTETED 

PARAMETER (1%=180,kx=60) 

DIMENSION S0(1X,4),¥0(60,4),AQ(IX,4),XP(1),YP(1) 


wewno 
Wow 


g 


TRIS SUBROUTINE UNWRAPS THE AIRFOIL AND STORES THE UNWRAPPED 
SURFACE IN ARRAYS AO AND SO. IT ALSO DETERMINES THE STRETCBING 
IN THE ETA DIRECTION. 


ZES88888 


g 


IMID = (11 +1) / 2 
DY = .8 / (33 - 2) 

D1lJe*2, 33 

Y © FLOAT(J-2) © DY 

¥O(J,1) © 1.25 * ¥ 7 (1. -¥ * ¥) 


1 


U = XIP{1) - XSING 
P(1) - YSING 


wn OO Ee. 


Rll © XP(I) ~ XSING 

Y11 = ¥P(1) - YSING 

ANGL = ANGL + ATAN2((U*Y21-Vex11), (U*X11+V*¥11)) 
R= SQRT(X11**2 + Yl1le#2) 

uo = XL 

vooYll 

R = SORT(R) 

AQ{I,1) = R * COS({.5 * ANGL) 

SO(I,1) = R * SIN(.5S * ANGL} 

WRITE (6,1000) ; 

WRITE (6,2000) (1,A0(I,1),SO(I,2),1 = 1 , TI) 
RETURN 

FORMAT(1X, ‘UNWRAPPED COORDINATES IN THE TRANSFORMED PLANE’) 
FORMAT(IS , 2620.8) 

END 
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SOURCE PROGRAM 


nse2d.f 


SUBROUTINE TABINT(XP,YP,XSING.YSING, 
Crocoussnssecern SRURKASSCHEERRA TS TEN ERH SCTE ERESSCERERKCESE 
ce * 
co SUBROUTINE TABINT * 
ee * 
CeeRnncceaanreeverarteeeeeenETee see EEO CORE RRAOSEL ENA SHOREEREMSOEEEEE 

PARAMETER (IX=160,kx=60) ; 

DIMENSION XP(1X),YP{IR),SO(IX) ,AO(IX) 

C = UP(1)} - XSING 

v= YP(1) - YSING 





vel. 

v=o. 

ANGL = 6. * ATAN(1.) 

po 11 =1,N ‘ 


X11 = XP(1) - ASING 

Y11 = Y¥P(I) — YSING 

ANGL = ANGL + ATAN2( (U*YL1-V*X11),(U*xXl14V*Y11)) 
R= SORT(XLL**2 + Y11 ** 2) 


R= SORT(R) 
AQ(I) = R * COS(ANGL * .5) 
1 S0(1) = R* SIN(ANGL * .5) 





61 
2 
3 


ry DX =(A0(N)-A0(1) 9/96. 

$862 | AOST = AQ(1) 

1863. po312*1, 97 

£864 XX = FLOAT(I-1) * DX + AOST 

“1865 | CALL TAINT(AO,SO,X%X,¥Y,N,3,NER,MON) 
1866: ZP(I) = XX * XR - YY * YY + XSING 
8671 93) YP(T) = 2. * 3X * Y¥ + YSING 

71868 | RETURN 

1869; END 


IK 
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SOURCE PROGRAM DATE ' 1/20/90 
nse2d.f TIME 49:46:04 am 


a tae ae ew. 












c AUAKHOCHERREDS TORRE KASTEN KEROEE 
ec * 
ce SUBROUTINE TAINT * 





ont 


zeeee 
wonauve 


c * 
CewrAaeeese#a4xwaeeese#enHeeSENRNSSCLERRATSOOETEUASO THE HEAEC SEER TASTE ENED 
PARAMETER (1X=180,kx7 50) 
DIMENSION 2TAB(1),FTAB(1),T(10) -C(10) 


MASA - ANES SUBROUTINE YOR POLYNOMIAL INTERPOLATION 
OF A TABULATED FUNCTION 


{3 


aang 


1 


IF(N-K) 12,1, 2 
NER = 2 
RETURN 
IF(K-$) 


a 
eeee ane 


SOaNAUSUN - 





@euwn Ww 
ia] 
" 
~ 

~ ew 


Se 


poe rei, NM 
IF(RTAB(I) - XTAB(I+1)) 9,11,10 
11 NER = 3 
RETURN 
3 = 3-1 
co T 8 
10 3 = J+1 
8 CONTINUE 
MON © 1 
IF(J) 12,6, 6 
12 MON = 2 
?7po13l*1,N 
IF(X - RTAB(I)) 14,14,13 
wage 
GO TO 18 
13 CONTINUE 

GO TO 15 
6po1l6r3*#1,N 
IF(2~XTAB(I)) 16.17,17 
g=y 
GO To 18 
16 CONTINUE 
iSJ=N 
183 = J - (K+1) / 2 
IF(J) 19,.9,20 
gel 
20M" 948K 
IF(M - N) 24,21,22 
2g=39-1 
GO TO 20 
2. KPl = K+ 1 
JSAVE = J 
DO 23 L = 1, KP1 
C(L) = X - XTAB(I) 
T(L) = FTAB(J) 
22.2 = Jel 
DO 2437 1.K 
Ie J+. 
25 T(E) = (CCT) PTL )-C(T) * TS) 7 (C(I )-C(T)) 
1211 
IF(I-KP1) 25,25,24 
24 CONTINUE 
FX = T(RP1) 
NER = 1 
RETURN 
END 


asin! 
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SOURCE PROCRAM 


nse2d.f 


SUBROUTINE SING(N2,N,X,2,2LE,YLE, TEA, TES, XSING , YSING, CED) 

f. COSCESAARHSSSHR ALE KTOSHETRR SRSRETASE STE RARASSHRERALHHSSHERESCRCSETNRAASET = 
jes : 

jer SUBROUTINE SING arse 
| ce : = E . 
Cate seceea nna se Saar Kss eee aR Nee RREASCOSEEEs UETKEV Se TIReAeeeeaene 


PARAMETER (I%=180,kx=60) 









c 
c 
c BIS SUBROUTINE COMPUTES SINGULAR POINT LOCATIONS. 
c 
DIMENSION X(2) , 2(2) 
NU = N2 
Nl = N24 1 ° 
N32 N2- 1 
X20 = 3(N1) 
Z1o= 2(K1) 
X20 © X(N2) 
220 = B(X?) 
X30 © (83) . 
23 = 2(K3) 
Dl = 42 +» 2-21 #* 2 
D2 © 22 ** 2-21 ** 2 
D3 = x2 - x1 
D4 = 22 - 21 
DS = X3 ** 2- x1 *# 2 
D6 = 23 ++ 2-21 e* 2 
D? = x3 - x2 
D8 = 23 - 22 


B= (D7 * ( Dl + D2) ~- D3*(DS+D6) )/(2.*(D7*D4-D3*D8) ) 

IF(ABS{D3).LT.ABS(D7)} GO TO 10 

A (D1 + D2 ~- 2. * B* DS) / (2. * D3) 

GO TO 20 
20 A = (DS + D6 - 2. * B* D8) / ( 2. * D7) 
20 CONTINUE 

R = SORT((X2-A)** 2 + (Z2-B)**2) 

AXLE = X(N!) 

YLE = Z(t) 

CED = X(1) - X(NU) 

AZ = (2(2)-2(2)) /(4(2) -— XQ1)) 

AL = (2(N)-2(N-1))/(X(N)-X(N-1)) 

TES = .5 * (AL + A2) 

TEA = A2- Al 

TEA = TEA « 57.29578 

XSING = (7-XLE}) /2. 

YSING = (B-YLE) / 2. 

RETURN 

END 








SOURCE PROGRAM DATE 1/20/90 7 | 
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BoP cee on: 


SOURCE TEXT 


SUBROUTINE AIRFOL(IT, IE, ILE) 
i, CASSHERATEASSTEAH LSAT HEERAKSHSEHRTSSSSEERATO SSH REKESHERERESSCSHEEESE OEE 
6tee. SUBROUTINE AIRFOL * 
Cet tees arere era case esssaheseePaae se ese Hugg st Hoansearsesentateteserette 
PARAMETER (1X%=180,kx=60} 
COMMON/GRID1/% (IX, KX) ,2(1X, KX) 
COMMON/DGRID/DT, IMAX, KMAX, ITEL, ITEU 
COMMON /Y SYM/I SYM 
DIMENSION S0(IX,4),AO(IZ,4),Y0(60,4),XP(IX) ,YP(IX), 
1E(IX). F( 1X), BO(49) 
DIMENSION XL(IX), RU(IZ), YL(IX), YUCIX), 
2 : RU(IX), YY(IX) 


DATA (BO(I),T=1,32)/2.,1.0414,1.0836,1.1270,1.1715,1.2175,1.2651, 
11.3145,1.3659,1.4199,1.4755,1.5349,1.5973,1.6636,1.7342,1.8099, 
21.8914,1.9799,2.0764,2.1829,2.3012,2.4341,2.5653,2.7597,2.9646, 
33.2206,3.5141,3.9019,4.4219,5.1687,6.3632,6.6809/ 


READ(S,*) IBMAX,ISYM 

IE( ISYM .EQ. 0 ) THEN 

DO 101 1=1,TBMAX 

READ(S,*) ZU(I), YD(I), YL(I) 
101 CONTINUE 

ELSE 

DO 102 1 = 1,IBMAX 
202 READ(5,*) XU(I), YU(I) 

po 103. I = 1,ITBMAX 
103 YL(I) = ~ YU(I) 

ENDIF 


feat 


DO 1000 I=1,1BMAX 
IU = I + IBMAX 
AX(IU) = XU(T) 
y(1u) = yu(1) 
IL = IBMAX - I+ 1 
AR1) © XUCTL) 
YY(1) = YL(IL) 
1000 CONTINUE 
c 


7 


vv 


Dt te neernes meee 


IBMAX2 = 2° I BMAX 
DO 1010 T=1,IBMAX2 
XP(I) * XX(1) 
YP(I) = YY(T) 
1010 CONTINUE 
FNU = IBMAX 
FNL = IBMAX 


CUCM 
VRWN = OO ONO RUN — O06 S) 


irvinyi 


THIS SUBROUTINE GENERATES SHEAR PARABOLIC C-GRID 

THE FOLLOSING SUBROUTINES ARE RELATED TO THE GRID GENERATION 
WRAP SING 
TABINT CLUSTR 
TAINT STRTC 


do g8I=1 , 32 
8 AOQ(I,1) = BO(I) 
II = IMAX 
JJ = KMAX+1 
we 31 
IE = 127 
TmPl «= Ir +21 
IIM1 = II - 2 
13 = I1 * 32 
TIJJ2 = 11 * ( JJ-2) 
ILE * (17 + IE ) / 2 
PI = 4. * ATAN(1.) 
NU = FNU 
NL = FNL 
N o= NU + NL 
SCALE = 1./ ( XP(1) ~ XP(NL) ) 
DO 3333 1=1,N 
XP(1) = XP(I) * SCALE 
YP(I) = ¥P(I) * SCALE 
3333 CONTINUE 
CALL SING(NU,N,2P,YP,XLE,ZLE, TEA, TES ,XSING , YSING, CHD) 
CALL TABINT(XP,YP,XSING,YSING,N) 
NBODY = I1£ + 1 - IT 
DO 6791 J = 1 , NBODY 
Lert-1 
E(IT+L) = XP(T) 
6791 F(IT+L) = YP(1) 
IEPL «= IE +1 
SLOPT = TES * .75 
po 438 I = IEP] , II 
Tl = 141 - ‘IE 
E(1) = AO(I1,1) 
OYI = 1. / 48. 
D #4. / 3. * (EI) - .25) 
F(1) = FCIE) + SLOPT * ALOG(D) / D 
L= pl - 1 
E(u) = Eq1) 
438 F(L) = F(IT) + SLOPT * ALOG(D)/D 

c WRITE (6,439) 

C439 FORMAT(2X,3H 1,19%,182,19X,1HY) 

c WRITE (6,37) (1,E(1),F(1),1 = 2, 11) 
CALL WRAP(11,JJ,XSING,YSING,E,F,S0,A0,¥0} 
pboid0sge2, JI) 
po1l0r=#1, 11 
R(I,J-1) = AOCT,2}*#2 - (SO(1,1)+¥0(T,1))**2 

10 Z(1,5-1) = 2. * AO(I,1) * (S0(7,1)+¥0(I,1)) 
RETURN 

37 FORMAT(1I5,2F20.8) 
END 
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SUBROUTINE CLUSTR(DS} 


ce SHARRARSSSTUERAS TREASMSHLETAASSSATTTCHSSSORLTAHSSSSRKESSSSERRARAEOSD 
ce  * 
ce SUBROUTINE CLUSTR . * 
ce * 
CeURAeS Ses aRKS See He AAHSESHTRRKSHSHHETKKSORHERHSSSEHEELRESSTERATHSTATETR 
PARAMETER (1X=180,kx=60) 

COMMON /GRID1/X (IX, RX), 2(IX, RX) 

COMMON /DGRID/DT , IMAX , KMAX , ITEL, ITEU 

DIMENSION S(60),%P(60},YP(60),R(60) 






FAIS SUBROUTINE CLUSTERS A GIVEN X,2% GRID SUCH TRAT THE FIRST POINT IS AT 


po100 1 = 1, IMAX F 
Sql) = 0. 
RP(L) = X(I,1) 
YP(1) = 2(1,1) 
pol0 k= 2), KAX 
EP(K) = X(1,K) 
YP(R) = 2(1,K) 
1o S(K) = SORT((XP(K)-XP(K-1))**2+ (YP(K)-YP(K-2))**2) 
14S(K-1) 
i ‘SUMDX = S(KMAX) 
| CALL STRTCH(SUMDX,DS,F1,KMAX, FACTOR) 
C | WRITE (6,200) 1, FACTOR 
RQ) = 0. 
DR = Ds 
DO 20 K = 2, KMAX 
RK) = R(K-1) + DR 
DR = DR * . ACTOR 
20 CONTINUE 
RLAST = 1. / R(KMAX) 
pO 30 KR = 2, KMAX 
Rl = R(K) * RLAST * S(RMAX) 
CALL TAINT(S,XP,R1,XP1 , KMAX, 3,NER, MON) 
X(1,K) = XP1 
CALL TAINT(S,YP,R1,YP2,KMAX,3,NER, MON) 
2(1,R) = YPL 
} 30 CONTINUE 


ana 


| 


rrr} 


— 
‘i 


BNAVEWN—-CVDNAVAWN-—OO BNO 


i 


NNN Ne 


we eee a a ee wr oe i we oe i 


100 CONTINUE 
ic WRITE (6,115) 
DO 110 I = 1. IMAX 


MNSNINININ SN 


: DX © X(1,2) - X(1,1) 

: DY = 2(1,2) - 2(1,2) 

; DN = SORT(DX*DX+DY*DY) 

ic WRITE(6,120) I , DX , DY , DN 

110 CONTINUE 

} RETURN 

125 FORMAT(5X,6HNORMAL,1X,8HDISTANCE,3B AT,4H THE,5R WALL,/ 
1,58 1,8X,2HDX, 8X, 2HDY, 2X, 2HDN,//) 

120 FORMAT(I5,3F10.5) 

200 FORMAT(I5,F10.5) 
END 


w 


datas 
1 BR I oe ore 
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SOURCE TEXT 
SUBROUTINE STRTCH(SUMDI,DX1,F1,N1,R) 


CACeceeeaAeoeaeaHeCHAEERRSAeeaKAEHSHEKERAKSAELERAHHEHETERASORARAAESSS 

* ¢ 

SUBROUTINE STRICH * 

* 

SECRET AHS SHERAASASSHTTLASSERERHAASORERASSSHRERASSHSUEEN ASCE TERE 


PARAMETER (1X*180,kx=60) 


THIS SUBROUTINE DETERMINES A GEOMETRIC 

PROGRESSION OF GRID SPACING BETWEEN 1 AND Nl SUH THAT 
BUNSDZ} EQUALS SUMDI. THB BATIO BETWEEN SUCCESSIVE 
SPACINGS I5 BR, 


« 
102 = 1, 50 
(R-1) * SUMDX - DxX1*(R**N-1) 
FP = SUMDX - DX1 * FLOAT(N) * R ** (N-1) 
RITER = P - F/ FP 
3P(1.E~-02.L1. RITER.AND.RITER.LT.1.) RITER = 1. 
T¥(1. .LT.RITER.AND.RITER.LT.100.) RITER=.01 
TP(ABS (R-RITER).LT. R*E1) GO TO 1 


fi 


sr val oh alah ahvel ahah arf ot ot oe ated 


3 


41 = DETOT/FLOAT(N1~1) 
RETURN 
Re RITER 
RETURN 
END 
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SUBROUTINE EDDY (CMUL) 
CRCCCERRAAACORNEH SSS ETERAAASEESEAGSSCRERTESO OU ERRESOOENRESOOCERERSS OE 
2 ice : : = 
jlo SUBROUTINE EDDY Eee © . 
ce oo * 
CHRTAMSeeOA RAKES ELE MRSA SeRRARA SSH HTRHASCOPRRAASSCOPRERSCRSEUNETEOSEERE 
PARAMETER (1X=180,kx=60) : 
COMMON/FLOW/Q1 (1X, KZ) ,02 (IX, KX) ,Q3( IR, KX) .Q4( IX, KX) 
COMMON /MUTUR/CMU (IX, KX) 
COMMON/SKINCF/CF(I2} 
COMMON/DGRID/DT , IMAX, KMAX, ITEL , ITEU 
COMMON /PAR/GAMMA , REYREF , ALFA, ALFA , REDFRE , AMINF,ALFAI 
COMMON/GRID1/2 (1%, KX). 2(1%, KX) 
COMMON/SPEED/ U(IX, KX), V( IX, KX) ,AA( IZ, FX) 
DIMENSION TIN(KX), POUT(KX) ,Y( KX), TRANS (IX), S(IX),UO( KX) ,VE(IX) 


CMUPP = 1000. * CMUL 

REDGE = KMAX 

ILE = ( ITEL + ITEU ) / 2 

CRORD = X(ITEV,2) ~ X(ILE,1) 

pO 100 IT = 2, IMAX - 1 

TPq{ ABB(X(i,1)).LT.{ ABS(X(ILE,1)) + 0.05 * CHORD 3 ) GO TO 100 
UDIF = 0. 

UMAX = 0. 

UMIN = 9999. 

MAX = 0.1E-06 

YMAZ = .1E-06 

PYMAX = 0.1E-06 

Yql) = 0. 

COMPUTE TAU AT THE WALL 

VET = U(1.2) - U(I,1) 

VET = V(I.2) - V(I,2) 

RXI = x(I-1,1) - X(1-1,1) 

OXI = BCi+1,1) - 2(1-1,21) 

XET > xCI,2) - 3. * ¥¢1,1) - 2(7,3) 

* 2(1,2) - 3. * 2¢1,1) - 2(1,3) 
* XXII 

+ Ux 


Phat ht Sr cl a aa nee at 
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. / (RXI * ZET - 2X1 * XET) 
(UET * XX1 + VET * 2XI ) * YAC 
= ANINF «© OMEGA / REYREF 
CF(1) = 2. * TWALL / (AMINF**2) 
FACT = SQRT(Q1(1,1) * ABS{TWALL)) *REYREF/ (26. *AMINF) 
DO 10 K = 2, KEDGE-1 
UXI = U(I-1,K) - U(I-1,R) 
V(I+1,R) - V(I-1,K) 
U(I. R41) - U(I,R-1) 
VqI,K+1) - V¢I,K-21) 
X(I-1,K) - X(T-2,8) 
2(1-1,R) - 2(I-1,R) 
A(I.R4L) = X(I,R-1) 
2(1,K+1) - 2(1,R-1) 
. 7 (MXIT * ZET = 2X1 * XBT) 
ABS (UET*XX1+VET*ZX1-UXI*XET-VXI*ZET) * YAC 
SORT(U(1,K)**2 + V(I,R}**2) 
AMAX1 (UTOT , UMAX ) 
AMINL (UTOT,UMIN) 
Y(R) = SQRT( (X(T, K)-(1,K-2))*924(2(1,R)-2(1,R-1) )**2)4¥(K-1) 
F = Y(K) > OMEGA 
IF((¥(R)*FACT).GT.20.) GO TO 31 
IF(1.GT.ITEL.AND.I.LT.ITEV) F = F * (1. - EXP(-¥(K)*FACT)} 
CONTINUE 


rag? pee ee i 
NRANUNN NNN 


RON NI IRN RR RI RI 


1F(F.GT.FYAX) TBEN 
FMAX = F 

YMAX = ¥(K) 
ENDIF 

FCT = Y(K) * FACT 

IF(ECT.GT. 20.) FCT = 20 

FCT = ABS( FCT) 
EL = .4° ¥(R) * (1. - EXP(-FCT)) 
TIN(K) = QL(I,R) * EL * EL * OMEGA 
TIN(K) © ABS(TIN(K)) 
CONTINUE 

UDIF = ABS(UMAX-UMIN) 

KSWTCE = 0 

FWARE = YMAX * FMAX 

Fl * 0.25 * YMAX * UDIF **2 / FMAX 
IF(FL.LT.PWARE) FWARE = £1 

Do 20 K = 2, KEDGE - 1 

FKLEB = 0. 

IF(ABS(Y(K)/YMAX) .LT.1.E+04) THEN 
FRLEB = 1. / (1. + 5.5 * (0.3 * Y¥(K)/YMAX) ** 6) 
END IF 

TOUT(K) = .0168 * 1.6 * Q1(I,K) * FWAKE * FRLEB 
TOUT(K) = ABS(TOUT(K)) 

IF(RSWTCH NE.0) GO TO 20 
IF(TIN(K) GT. TOUT(K)) KSWICB © K - 1 
CONTINUE 

DO 30 K © 2, KEDGE - 1 
IF(K.LE.RSWICE) THEN 

CMU(1,&) = ABS(TIN(K)) 

ELSE 

CMU(I,R) = ABS(TOUT(K) ) 

END IF 

CONTINUE 
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FROZE TBE VALUE OF EDDY VISCOSITY TO AN UPPER LIMIT 


Rel 

Reka 

If(K.GT.REDGE) GO TO 736 
IF(CMU(1,£).LE.CHUPP) GO TO 734 
CMU{1,%) = CMUPP 
grersi 
If(K.LE KEDGE) GO To 735 

CONTINUE 

CONTINUE 

RETURN 

END 
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SUBROUTINE RESI 


CASSHHERRASSHRAREASSHARTERSSOCERRAESSOEREREHSETURECHETENEAHA CIEE RAATES 
. 
SUBROUTINE RESI * 
t 
CeRMAHH SHAR AAHSSHRARKSSCSHERERAHSSHERARAACOHERAALA SHA LERHGS TER TTAHO TEER 
PARAMETER (I%=160,Kx=60) 
COMMON /FIX/OMEGA , EDOT 
COMMON/PERTR/DQ1( IK, KZ) ,DOQ2(IX, EX), DQ3(IX, KX) ,DO¢ (IX, KX) 
COMMON/GRID1/X (IX, KX), Z(1X, KX) 
COMNON/DGPID/DT, IMAX, KMAX , ITEL, ITEU 
COMMON /FLOW/Q1 (1%, KX) ,Q2( IX, KX) ,Q3( IX, KX) ,Q4( IX, KX} 
COMMON/SPEED/ U(IX,KX),V(IX,KX), AA(IZ, KX) 
COMMON /PAR/GAMMA , REYREF , ALFA, ALFA1 , REDFRE , AMINF , ALFAI 
COMMON/COFRES/ RHS(IX,4) 
ATAV(I,K) = OMEGA * 2(1,K) 
YTAU(I,K) = - OMEGA * X(I,K)- MDOT 
EIS SUBROUTINE COMPUTES THE RESIDUAL ON THE RIGHT HAND 
BIDE ARISING PROM THE EULER- PART OF THE GOVERNING EQUATIONS 


ELUX GERMS WITHIN THE II- DERIVATIVE 
DO 100 K = 2, KMAX- 1 
po 10 I= 1, IMAX 
UCON = UCI,K) © (2(1,K+1)-2(1,K-1)) 
Li- VEILR) * (XC1,R+1)-2(1,K-1)) 
UCON « 0.25 * DT * UCON 
RIT © - ATAUCI,K) *(2(1,K+1)-Z(3,K-1)) 
1 + YTAU(I,R) * (X¢1,K41) - X(1,R-1)) 
RIT = RIT * DT * 0.25 
UCON = UCON + XIT 
RBS(I,1) = UCON * Q1(I,K) 
P= (GAMMA-1.) * (Q4(I,R) - .5*QU(Z,K)*(U(I,K)**2 + V(I,K)**2) ) 
RAS(I,2) © Q2(1,K) * UCON + P * DT * 0.25 * (Z(I,K+l) - 2¢1,8-1)) 
RHS(T,3) = Q3(1,K) * UCON - P * DT * 0.25 * (X(1,K+1)-X(I,R-1)) 
RHS(I,4) = UCON * (Q4(1,K)+P) - X1T * P 
CONTINUE 
DO 1l I= 2, IMAX~1 
DO1{1,K) = DO1(I,K) - RBS(I+1,1) + RHS(I-1,1) 
DQ2(1,K) = DO2(I,R) - RBS(I+1,2) + RBS{1-1,2) 


: pa | { 


DO3(1,K) = DQ3(1,K) - RBS(1+1,3) + RBS(I-1,3) 
© DO4(I,R) - RES(I+1,4} + RHS(I-1,4} 


peryporrerenrererereaer ra ad 
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DOQ4S(I,8) 
CONTINUE 


FLUX TERMS WITBIN THE ETA- DERIVATIVE 


DO 200 I = 2, EIMAX - 1 

DO 20 K = 1, KMAX 

VCON = U(I,K} * (2(3-1,8)-2(I41,K)) 
1 4V(I,K) * (X(1+1,R)-X(1-1,K)) 

VCON = VCON * 0.25 * DT 

ETAT = -XTAU(I,R) * (2(1-1,K)-Z2(141,K)) ~ YTAU(I,R)* 
1 (RC141,R)-X(1-1,K)) 
ETAT © ETAT * 0.25 * DT 

VCON = VCON + ETAT 

RHS(K,1) * VCON * Q1(I,R) 

P = (GAMMA-1.) * (Q4(1,K) - 0.5 © QUCTLR)*(U(T,R)e*2 + V(I,R)**2)) 
RAS(K,2) = VCON * Q2(1,K) +P * DT * .25 * (2(3-1,K)-Z(141,K)) 
RAS(K,3) » VCON © Q3(1,K) + P * DT * .25 * (X(14+2,K) - X(I-1,K)) 
RHS(K,4) = VCON * (Q4(1,K)+P) - ETAT * P 

CONTINUE 

DO 21 R= 2, KMAX ~ 1 

DQL(I.R) = DQL{I,K) - RBS(K+1,1) + RES(K-1,1) 

pDQ2(1,R) DQ2(1,K) - RBS(K+1,2) RBS(K-2,2) 

DO3(1,R) = DQ3(I,F) - RBS(K+1,3) * RHS(K-1,3) 

DQ4(1,R) = DO4(3,K) - RHS(K+1,4) + RHS(K-1,4) 

CONTINUE 

FORMAT(216 ,4£14.6) 

RETURN 

END 
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SUBROUTINE ROTGRID( IMAX, KMAX, DALFA) 

Crceuugsrreecoewussaressenns SRAETAKKHSSEARARSESSSHETATSE SE SREKKROSCSERRACESSE 

e 
SUBROUTINE ROTGID e 

e 

Ceeeeesceeweer awe cee el Asse eH #&RRHHSHTERAHSSSERATELACTOLERASHHEEERESS CRETE 

PARAMETER (1X=180,kx=60) 

ROTATE GPID IN THE CLOCKWISE DIRECTION BY AN AMOUNT DALFA 

COMMON /GRID1/X(1X,KR)},2( IX, RX) 

CA = COS(DALFA) 

SA = - SIN(DALFA) 

DO 20 K= 1 , KMAX 

DO 20 1 = 1, IMAX 

X1 = X(1,K) 

Zi = 2(1,%) 

X(1,K) = 21 * CA - 22 © SA 

Z(I,R) = Zl * CA 4+ XL * SA 

RETURN 

END 
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SUBROUTINE CPPLOT 


375 } COMMA Pe eee REAReeeeHERFECERERAAOSERETHTOSOLURISOOSERREL ESOS ETRISOSORES 
Le PARAMETER (119180, kx=60) 

78 ‘RIS SUBROUTINE PLOTS CP AT EQUAL INTERVALS IN THE MAPPED PLANE 
COMMON/SURF/PSUR (IX) 

DIMENSION KODE(4), LINE(90),X(IX),¥ (IX) 
DATA KODE/1H ,1H+,181,18*/ 
WRITE ( 6 . 2) 

2 PORMAT(SOROPLOT OF CP AT EQUAL INTERVALS IN THE MAPPED PLANE/ 
2 1080 zx , 108 xvc. 108 CPL ,108 CPL 
CPO = (1. + .2 * FMACH **2) ** 3.5 - 1. 

CPO = CPO / ( .7 * FMACH **2) 
KO = 30. * cpO + 4.5 

IMIN = (32-1%1)/2 + I1 

ILOW = 2:* IMIN 

CBD=X(11) - X(IMIN) 

polz1= 1, 90 

LINE(I) = KODE(1) 

DO 34 1 = IMIN , I2 

KR = 30. * (CPO - PSUR(I}) + €.5 
KL = 30. * (CPO -PSUR(ILOW-I)) + 2 5 
IF(K.LT.1) R= 4 

IF(KI.1T.1) Kl = 1 

IF(K.GT.90) K = 90 

IF(K1 .GT. 90) Kl = 90 

LINE( KO) = KODE(3) 

LINE(K) = KODE(2) 

LINE(K1) = KODE(4) 

MOC = (X(1) - X(7MIN)) / CBD 
WRITE (6,610) X(1),XOC,PSUR(ILOW-1),PSUR(1),LINE 
LINE(K1) © KODE(1) 

LINE(K} = FKODE(1) 

RETURN 

FORMAT(4F10.4) 
FORMAT(4F1C.4,90A1) 

END 
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nse2d.f 
SOURCE TEXT 
2412: SUBROUTINE OUTRVC(REREAL) . 
2413 } CeeeeaaateecekansaeeeURAHAOKSaHRAASSCCEREMOSSERARERS HERR ESS ON SERRE S ES 
"2414 } ce * 
_24is = SUBROUTINE OUTRVC . 
2416 | ce * 
7 | CemanscorverreceowranesereneeneCenenheeeeseREteeeessanseersnaneeensEn 
PARAMETER (1X=180,kx*60) 
COMMON /PLOT/TITLE(10) ,NSTPT, RESD( 3000) ,RES ,CLB( 3000) .CDPH( 3000) 
COMMON /DGRID/DT, IMAX, KMAX , I TEL, ITEU 
COMMON /FLOW/Q1 (IX, KX) ,02(1X, KX) ,03 (1X, ER) ,Q4( 1X, KX) 
COMMON /PAR/GAMMA ,REYREF , ALFA, ALFA], REDFRE , AMINE , ALFAI 
COMMON /SURF/PSUR(IX) 
COMMON /SKINCF/CF(IX) 
COMMON /MUTUR/CMU (IX, KX) * 
COMMON/TITL/ITITLE 
CHARACTER ITITLE*80 
2428. DIMENS ON FY(IZ,K2) 
a8 DUM = 0. 
CODE = 0. 
431 IRES = IFIX(RES) : 
432 | ALFAD = ALFA * 45. / ATAN(1.) 
2433 CHUL = AMINF / REYREF 
2434 | pO 10 Ke 1, KMAX 
2435 | po 101 = 1, IMAZ 
2436 | FY(I,K) = 0. 
2437 | 10 CMU(I,R) © CMU(1,E) / CMUL 
2438: WRITE(4) ITITLE 
2439 | WRITE(4) IMAX,KMAX, TEL, ITEU, AMINF, ALFAD, REREAL, DUM,NSTPT, GAMMA, 
"2440 | + CODE, RES. DUM 
244) | WRITE(4) Q1,07,03,04 
2442 | WRITE(4) RESD 
2443 WkKITE(4) PSUR 
2444 WRITE(4) CF 
| 2445 | WRITE(4) CLH,CDPR 
2446! WRITE(4) CMU 
2447, WRITE(4} FY 
2448 WRITE(4) TK 
2449 | WRITE(4) TE 
2450 RETURN 


“2451! END 
| 
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SUBROUTINE OUTNT(NOUT,NSTPP, TIME , REREAL) 
PARAMETER (IX=160,kx~60} 

COMMON/GRID1/X{1X,KX),2( 1%, KX) 
COMMON /PAR/GAMMA , REYREF , ALFA, ALFA] , REDFRE , AMINF , ALFAI 
COMMON/DGRID/DT, IMAX , KMAX , ITEL, 1 TEU 


COMMON / FLOW/Q1 (IX, KX) ,Q2(12, KX) ,Q3(IX, EX) ,Q4(1E, KX) 
ITN * NSTPP 


PIs 


4.*ATAN(1.) 


ALPAD*ALFA* (180. /PI} 


REWIND 
WRITE 


IF( IT .£Q. L*NOUT ) THEN 
31 
(31) IMAX , RMAZ 


WRITE (31) ( ( X(I,K}, I=1,XMAX ), K=l, MAX 
WRITE (31) ( ( Z(1,K), I=l,IMAX ), Kal, KMAX 
WRITE (31) IMAX , KMAX 

WRITE (31) AMINF, ALFAD, REREAL, TIME 


WRITE 
WRITE 
WRITE 
WRITE 
REWIND 
END 1F 


(31) ¢ ( Q1U(I,K), Iml,IMAX ), K=1,KMAX 
(31) ¢ ¢ Q2¢1,K), I=1,2IMAX ), Kel, KMAX 
(32) € (¢ Q3(1,K), Isl, IMAX ), K=1,KMAX 
(32) ¢ ¢ Q4(I,K), I=2,IMAX ), K-1,KMAX 
31 


IF( ITN .EQ. 2*NOUT ) THEN 


REWIND 32 

WRITE (32) IMAX , KMAX 

WRITE (32) ( ( X(I,K), I™1,IMAX }, Rei, KMAX 
WRITE (32) { ¢ 2(1,R), I*2,IMAX ), Kel, KMAX 


WRITE 
WRITE 


(32) IMAX , KMAX 
(32) AMINF, ALFAD, REREAL, TIME 


WRITE (32) ( ( Q1L(I,K), I=1,IMAX ), K=1,KMAX 
WRITE (32) { (¢ Q2(I,K}, T=l,IMAX ), R=l, KMAX 
WRITE (32) ( ( Q3(I,R), I=1,IMAX ), Rel, KMAX 
WRITE (32) ( ( Q4(I,K), I=1,IMAX ), Rul, KMAX 


REWIND 
END IF 


32 
IF( ITN .EQ. 3*NOUT ) THEN 


REWIND 33 

WRITE (33) IMAX , KMAX 

WRITE (33) ( ( X(1,K), I=1,IMAX ), Kel, KMAX 
WRITE (33) ( ( 2(1,K), T©1,IMAX ), Kel, KMAX 
WRITE (33) IMAX , KMAX . 
WRITE (33) AMINF, ALFAD, REREAL, TIME 

WRITE (33) ( ( QL(I,K), T=1,IMAX ), Rel, KMAX 
WRITE (33) ( ( Q2(I.K), I=1,IM'X ), Kel, KMAX 
WRITE (333 ( ( Q3(I,8), Isl,IMst ), K=1,KMAX 
WRITE (33) ( ( Q4(1,K), I=1,IMAX ), Kel,KMAX 
REWIND 33 

END IF 


IF({ ITN .EQ. 4*NOUT ) TEEN 


REWIND 34 

WRITE (34) IMAX , KMAX 

WRITE (34) ( ( X(1,K), I=1,IMAX }, K=1,KMAX 
WRITE (34) ( ( 2(I,F), I*1,IMAX }, K=1,KMAX 
WRITE (34) IMAX , KMAX 

WRITE (34) AMINF, ALFAD, REKEAL, TIME 

WRITE (34) ( ¢ QUEI,K), I=1,IMAX ), K=1,KMAX 
WRITE (34) ( ¢ Q2(1,%), 1#1,IMAX ), K=1,RMAX 
WRITE (34) ( ( Q3(1,F), T=1,IMAX ), K=1,KMAX 
WRITE (34) ( ( Q4(1,K), I=1,IMAX ), K=1,KMAX 
REWIND 34 

END IF 


IF( ITN .£Q. S*NOUT ) THEN 


REWIND 35 
WRITE (35) IMAX , KMAX 
WRITE (35) ( ( X(I1,K), I=1,IMAX ), K=1,KMAX 


WRITE (35) ( ( Z2¢1,K), T=1,IMAX ). K=1,KMAX 
WRITE (35) IMAX , RMAX 

WRITE (35) AMINF, ALFAD, REREAL. TIMF 

WRITE (35) ( ( QLCI,K). T#1,IMAX ), K=1,KMAX 


WRITE 
WRITE 
WRITE 
REWIND 
END IF 


REWIND 
WRITE 
WRITE 
WRITE 
WRITE 
WRITE 
WRITE 
WRITE 
WRITE 
WRITE 
REWIND 
END IF 


(35) ( ¢ Q2¢1,K), T=1.IMAX ), K=1, RMAX 
(35) ¢ ¢ Q3(2,.K), Tel,IMAX ), R=}, KMAX 
(35) (¢ ( QS(T,R), I*1,ITMAX ), R=1,KMAX 
35 


1F( ITN .2Q. 6*NOUT ) THEN 

36 

(36) IMAX , KMAX 

(36) ¢ ¢ XC1,K), Tel,IMAX ), K=1,KMAX 
(36) ( ( 2(1,R), Tel, IMAX ), h¥1,RMAX 
(36) IMAX , RMAX 

(36) AMINE, ALFAD, REREAL, TIME 

(36) ( ¢ QECI,K), T=l,IMAX }, R=1,KMAX 
(36) ( ( Q2(1,R), Tet, IMAX }, Kel, KMAX 
(36) ( ¢ Q3(1,K), T©1,IMAX ), K=1,KMAX 
(36) ¢ ¢ Q4(E,R), 1=#1,1MAX }, K=1,RMAX 
36 


IF( ITN .£Q. 7*NOUT ) THEN 


REWIND 37 

WRITE (37) IMAX , KMAX 

WRITE (37) ¢ ( ACI KR), Iel1,IMAX ), Ke1,KMAX 
WRITE (37) (¢ ( 2(1,K8), Iml,IMAX ), Rel, KMAX 
WRITE (37) IMAX , KMAX 

WRITE (37) AMINE, ALFAD, REREAL, TIME 

WRITE (37) ¢ ( QL(I,K), T#t,IMAX ), K=1,KMAX 


REWIND 
END IF 


REWING 
WRITE 
WRITE 
WRITE 
WRITE 
WRITE 
WRITE 


37) ¢ ¢ Q2¢3,F), Tal, IMAX ), K=1,KMAX 
637) ( ¢€ Q3CT. KR). Ts], TMAX ), Rel, EMAX 
(37) ( € Q4CT,R), Th, IMAX ), Rel, RMAX 
37 


IF( ITS .EQ 8*NOUT ) THEN 

38 

(38) IMAX , KMAX 

(38) ( ( ECI,K), Tel, IMAX ), Kel, RMAX 
(30) ¢ ( 21,8), Tel, IMAX ), Fel, KMAX 
(38, IMAX , RMAX 

(38) AMINF, ALFAD, REREAL, TIME 
(38) ¢ (OGTR), 181, 2MAX ), ReL, MAX 


WRITE (38) ¢ ¢ Q2(1,R). I#1,IMAX ), F=i. RMAX 
WRITE (38) ¢ ( Q3(3.F), I=1,IMAX ), F-1,RMAX 
WRITE ¢38) (¢ ¢€ Q4(I.F}, T#l,IMAX ), Fr1,EMAX 
REWIND JF 
END IF 

TF, ITS .EQ) 9*NOUT ) TREN 
REWIND 39 
WRITE (39, TMAX , RMAX 
Witte 639) ( ¢ XCT,F). T-1,IMAX ), Fel. FMaAY 
WeiTE (39) ¢ ¢ ZCILE), Tel, IMAX ), Kel, REMAX 
WRITE (39: IMAX . EMAX 
WRITE (39 AMINE. ALFAIS. REPEAL. TIME 





wee 


~~~ 


Perereres 
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== 




















2572 | WRITE (39) ( ( Q2 
i i (I,K), 1=1,IMAX 
ae pis (39) ¢ ¢ Q2(3,K), T=l,IMAX 4 } 
28. mize (39) € ¢ Q3(7, 8), T=1,IMAX ), ) ; 
aaah i wire 439) ( ( Q4(I,K), 1*1, IMAX ), K=1,KMAX ) : 
“2577 | IF 
: ( IT .£Q. 10*NOU 
2398 : REWIND 40 pane 
3378 WRITE (40) IMAX , KMAX 
F550: WALTE (40) ( ( X€1,K}, 1=1,IMAX ), R=l,KMAX } 
geeks WRITE £90) ( ( 2(2,R), I=1,IMAX ), Kel, KMAX 
2583 | WRITE pee a inday 
$383 MRITE (40) AMINF, ALFAD, REREAL, TIM® 
RBs Gite ce € € Q1lC1,K), I=1,IMAX ), F=1,KMAX ) 
3586 | WRITE (40) ( COStIR), IeLcIMAK yy RCLGEMAX 
‘ Ky), 1*1, )}, K=1,EMAX ‘ 
B88 WITE (40) ( ( Q4(1,K), IS1,IMAX ), Kel, KMAE : 
89 | If( ITN .£Q. 11*NOU 
$390 | REWIND 41 sa cas 
592 | weize (al) Ct XC ma 
392. ( X€1,K), 1*1,IMAX ), R= 
3383 | WRITE (41) (| 2(1-), TeLsIMAX " rot RAN : 
2595 | WRITE Oh) Mane. 
383 | WRITE (41) AMINF, ALFAD, REREAL, TIME 
Beeei WRITE a1) ( ( Q1(1,K), 1=1,IMAX ), K=1,KMAX ) 
S38 ERITE (41) ( ¢ Q2¢1,K), I=1,IMAX }, R=1,KMAX ) 
t99 | RITE ore ( ¢ Q3¢E,%), I=1,IMAX ), K=1,KMAX ) 
2600 | baits ee > © ¢ Q4(1,K), I*1,IMAX ), K¥1,EKMAX ) 
a : If( ITN .EQ. 12% 
3602 | ces r Q. 12*°NOUT ) THEN 
$603 | WRITE (42) IMAX , KMAX 
Sos write (42) ¢ ¢ RCI,K), 121,IMAX ), K=1,KMAX ) 
2603 were 2) ( ( Z(1,K), T=1,IMAX ), el,KMAX ) 
2607 | WRITE CP OnE, 
$607 | WRI (42) AMINE, ALFAD, REREAL, TIME 
g608 ; mRITE (42) € ( QLQI,K), T81,IMAX ), K=1,KMAX ) 
rite WRITE (42) (¢ ¢€ Q2¢1,K), I#1,IMAX ), K=1,KMAX ) 
£8l0: BRITE (42) ( ( Q3¢1,K), 1*1,IMAX ), K=1,KMAX )} 
38th. me E ) ( ¢ Q4¢I,K}, I=1,IMAX ), K=1,KMAX ) 
“2613; IF( 1TN *NoU y 
zea penate Py EQ. 13*NOUT ) THEN 
2813 WRITE (43) IMAX , RMAX 
2616 waite (43) ( ( 01,K), I=1,IMAX ), K=1,KMAX ) 
abr: WRITE (43) ( ( 2¢1.K), 1=1,1MAX }, K=1,RMAX ) 
"2619 | WRITE a ear Hanes RE! 
; “INF, , REREAL, TIME 
beer wit uy) ( ( Q1(1,K), 1=1,IMAX ), K=1,KMAX ) 
=peel WRITE (4 ) ( ¢ Q2¢1,8), T21,IMAX ), K=1,KMAX ) 
se eats 43) ( ( Q3(1,K), T=1,IMAX ), R=1,KMAX )} 
2624 iD IF ) € ¢ Q4CI,KR), T=1,IMAX ), K=1,KMAX ) 
: IF( ITN .EQ. 14*NOU : 
Bee Cocgrel Q. 14*NOUT ) THEN 
| £627: WRITE (44) IMAX , KMAX 
$629 | metre (44) ( ( XC1,K), T22,IMAX ), K=1,KMAX ) 
nee WRITE (a4) ‘ ( 2(1.R), 1=1,IMAX ), K=1,RMAX ) 
263) | Terre tat Mae - REA 
ee eerie (44) AMINF, ALFAD, REREAL, TIME 
ae es (44) ( ( QL(I,X), 1=1,JMAX ), K=1,KMAX  ) 
ere are ia) ( ¢ Q2(1,R), Isl, IMAX ), Rel, KMAX  ) 
aoe wate (44) ( ( Q3¢1,K), I=1,1MAX ), K=1,RMAX  ) 
5636 kits Ete) ( ( Q4(1,K), T=1,IMAX ), K=1,KMAX ) 
2637 | RETURN 


2638 | END 
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PROGRAM AIRFGRID 

PARAMETER (1%=300, KX=100) 
COMMON/GRID1/X(IX, KZ) ,2(14, KX) 
COMMON/DGRID/DT, IMAX, KMAX, I TEL, I TEU 
LOGICAL VISCOUS 
READ(5,10) 
READ(S,100) IMAX,KMAX 
READ(5,10) 
READ(5,100) ITEL,ITEU 
READ(5,10) 
READ({5,200) Viscous 
READ(5,10) 
READ(5,300) DMIN 
READ(5,10; 

READ(5,*) AORAT , AOEXP , sdispl 

ILE = ( ITEU + ITEL ) / 2 

uP = ( ITEU - ITEL ) / 2 

WRITE(6,1000) IMAX,KMAX,ITEL,ITEU,ILE,IUP,DMIN,AGRAT 
1000 PORMAT( ‘Imax = ',710,10x,’kmax = ‘',310,5x,/, 


1 ‘Trailing edge lower =',110,/, 
2 ‘Trailing edge upper =',110,/, 
2 ‘Leading edge =',110,/, 
2 ‘lupper = Ilower =',110,/, 
3 *Distanpse of the first point =',f£20.10,/, 
4a "Streching ratio =',£20.10,/) 


CALL AIRFOL( AORAT , AOEXP , sdispl } 
IF( VIscOUS ) CALL CLUSTR( DMIN ) 





WRITE(6,1100) 
1100 FORMAT(//,'GRID BOUNDARIES',/) 
RTEOUT = ABS( X(1,1) - X(ITEL,1) } 
RLEIN © ABS( X(ILE,KMAX) - X(ILE,1) ) 
i IUP = «ILE + (¢ ITEU - ILE ) / 2 
i 


eer 


MN, 


RUP = ABS( Z(IUP,KMAX) - 2(IUP,1) } 
WRITE(6,1290) RTEOUT , RLEIN , RUP 


i 1 79x, 'Distanse * leading edge and inflow =',f20.10,/ 
i 2 ,5x,'Distanse of the body from the upper boudary=',f20.10,/) 

39: REWIND 21 

_.40: WRITE (21) IMAX, KMAX 

we Ah | WRITE (21) (( X(I,K), I21,IMAX), K*l,KMAX ), 

wn. 42: 1 (€ 2¢1,8), I*1,IMAX), K=1,KMAX ) 

3 STOP 

as: 10 FORMAT(1X: 


45 100 FORMAT(215) 
46 200 FORMAT(3L3) 
47 300 FORMAT(4F.0.0) 
48 | END 





SOURCE TEXT 


1200 FORMAT(5X.'Didtanse between trailing edge ans outflow =',f20.16,/ 
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airfgr.f 


SUBROUTINE METRIC 
CeCe eex&K Hee Kl ek HHH URE HOSEHARKAASTUARALSSSTRAKRASSHERATHEVCSSREEKASOES 
ee * 
SUBROUTINE METRIC : * 
* 
RRA HSSSHSEAARHSSHECRASCSSLAARHHSSSSTHRAHSSSLARAKH SHE C KARATE HKRHASESTEREE 
PARAMETER (1X*300,KX=100) 
COMMON / FIX/OMEGA , EDOT 
COMMON/DGRID/DT, IMAX, KMAX, ITEL, ITEU 
COMMON/GRID1/2( IX, RX) .2(12,KX) 
COMMON/GRID/YACOB( IX, KX) 
COMMON/MTRIX/XIX (IX, KX), RIZ(IN, KX), ZETAN(IN, KX), ZETAZ(IX, KX), 
AXIT(IX, RX), ZETAT( IX, KX) 


VE i: 
Gey a tas 


ype 


; igielo'g 


TRE SUBROUTINE METRIC COMPUTES TEE METRICS IN ALL THE TWO DIRECTIONS AND 
THE UNSTEADY COEFFICIENTS ETAT ZTC. 


DO 1000 K = 1 , RMAX 

DO 1000 1 = 1 , IMAX 

XTAU = OMEGA * 2(1,K) 

YTAU = OMEGA * (-X(I,K))- BDOT 

PRESENT SET UP IS POR FLOW PAST AK AIRFOIL. 


IP(I.£Q.1.OR.1.£Q.IMAX) GO TO 10 

EXT = 650 = (X(141,K)-X(I-1,K)) 

2X1 © 150 * (2(141,K)-2(2-1,R}) 

GO To 15 

IF(I.EQ.IMAX) GO TO 16 

EXE = 1.0 © (X(2,R) - X(1,K)) 

2X1 2 1.0 * (2¢2,K) - 2(1,K)) 

GO TO 15 

XXI © 1.0 * (X(IMAX,K) - XCIMAX-1,K)) 

2X1 = 1.0 * (ZCIMAX,K) ~ 2(IMAX-1,K)) 

CONTINUE 

IF(K.EQ.1.OR.R.EQ.KMAX) GO TO 17 

AZET = .5 *(X(1,R+1)-X(1,K-1)) 

ZZET = .5 *(2(1,K+1)-2(1,R-1)) 

GO TO 20 

IF(K.EQ.RMAX) GO TO 18 

RET = 2. * X(1,2)-1.5 * K(1,1) ~ .5 * 2(1,3) 
2ZET = 2. * 2(1,2) ~ 1.5 * 21,1) - .S * 2(1,3) 
GO TO 20 : 
XZET = 1.5 * X(1,RMAX)-2.* X(1,RMAZ-1)+.5°X(1,KMAX-2) 
ZZET = 1.5 * 2(1,RMAX)-2.* 2(1,RMAX-1)+.5*Z(1, RMAX-2) 
CONTINUE 

YACOBI = XXI * 2ZET - XZET * 2X1 

YACOB(I,®} = 1. / YACOBI 

XIX(I,K) = 2ZET * YACOB(I,K) 

XIZ(1.R) = -XZET * YACOB(I,K) 

XTAU = OMEGA * 2(1,K) 

YTAU = - OMEGA * X(1,K)- BDOT 

MIT(I,K) = - XIX(T,K) © ATAU - AIZ(I,KR) * YTAU 
ZETAX(I,R) * -ZXI * YACOB(I,R) 

ZETAZ(1,K) = XXI * YACOB(1,K) 

ZETAT(1,K) = - ZETAM(I,K) * XTAU - ZETAZ(1,K) * YTAU 
CONTINUE 

RETURN 

END 


aeiat. 


lettadolai 


| 


cakes seven 


0010) cn! on! cn! o0! cn! 
N—OOONAUM 








SOURCE PROGRAM DATE 1/20/90 


airfgr.f TME 11:45:57 am 
ities i a SM a ates o 


SUBROUTINE WRAP(II,JJ,ISING ,YSING ,XP,YP,S0,A0,Y0) 


RUAOSCRERRAAHSORETAHSSSERTKRESTHORAKESSSHNTSESECCERRROTEE 


SUBROUTINE WRAP 


SARTCSOHAETRASCHERAASRESERRALSSCEEURAHEOSRRAERSSHTAEKESCAERARASTENTIE 


PARAMETER (IX©300,K1=100) 
DIMENSION S0(1X,4),YO(IX,4),AOQ(IX,4),XP(1),¥P(1) 


THAIS SUBROUTINE UNWRAPS THE AIRFOIL AND STORES THE UNWRAPPED 
SURFACE IN ARRAYS AG AND 50. IT ALSO DETERMINES THE STRETCHING 
_ IN THE ETA DIRECTION. 


IMID = (11 +1) / 2 

DY = .8 / (33 - 2) 

polg=2, 33 

Y = FLOAT(J-2) * DY 

yO(J,1) = 1.25 * ¥/ (1. -¥* ¥) 

Yoql , 1) = - ¥0(3,1) 

PI = 4. * ATAN ( 1.) 

ANGL © PI + PI 

U = XP(1) ~ XSING 

Vo= YP(L) ~ x5iNG 

vel. 

v=0 

TIM = Ir ~ 1 

po2r=1, 11 

X11 = XP(1) - XSING 

Y11 = YP{I) - YSING 

ANGL = ANGL + ATAN2((U*Y11-V*X11),(U*Z11+V*¥121)) 
R= SQRT(X11**2 + Y11+*2) 

uo xll 

vo s®Y1i 

R = SQRT(R) 

AO(I,1) = R * COS(.5 * ANGL) 

$0(1,1) = R * SIN(.5 * ANGL) 

WRITE (6,1000} 

WRITE (6,2000) (1,A0(1,2},S0(I,2},1 ~ 2 , 11) 
RETURN 

FORMAT(1X, ‘UNWRAPPED COORDINATES IN THE TRANSFORMED PLANE‘) 
FORMAT(IS , 2F20.8) 

END 
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airfgr.f 










B 7 i 





SUBROUTINE TABINT(XP,YP,XSING, YSING ,N, sdisp1) 
CACSCHAAAACSSLER AROSE SET AERSSHREAATOTSRETRASHSHATAROCHESRARG SE LELTES ES 
e 
SUBROUTINE TABINT * 
* 
CeRanKeSevene Kae eee anASOeaRRKSS SHU KRKAESEHARAHOSSAREEKAS ER ATNHSOTERRE 
PARAMETER (IX=300,K%=100) 
DIMENSION XP(IX), YP(IX) ,SO(IX) ,AO( IX) ,DTHXL (IX) ,DXW( IX) ,DX(IX) 
COMMON /DGRID/DT, IMAL, KMAX, ITEL, ITEU, ILE 
PI = 4.*ATAN(1.) 
U = XP(1) - XSING 
¥v = YP(1) - YSING 


UTE 


aoe 


tales 
ALISE 


Qantun—oo 


O00" 


. * ATAN(L.) 
1,N 
Z1l = XP(1) - XSING 
YIL = YP(I) ~ YSING 
ANGL © ANGL + ATAN2((U*®Y11-V*xX11), (0*X11+V*¥11}) 
Roe SORT(R11**2 + Yli ** 2) 


ne) Oa 
SISEaF 


R= SORT(R) 
AO(I) = R * COS{ANGL * .5) 
SO(I) = R * SIN(ANGL * .5) 


Re 


“AOL = AQ(N) ~ A0(1) 
DIWs = 0.0 


AOST = A0(1) 
idiv = 1+ ( iteu - itel ) 


aah 


TAaDX «= 0. 
DTBXI(1) = 0. 
DTEXI(2) = pit SDISPL/2. 
DTEX = (1.~ SDISPL } * PI / (IDIV-1) 
DO 9 i= 3,IDIV 
DTBXI(1) = DTBXI(I-1) + DTEX 
DO 10 I = 1,1DIV 
THDX = (I~1)*DTEX 
if( 2 .eq. 2 .or. i .eq. idiv ) then 
if( 4 .eq. 2) thax = .5*sdispl * pi + thax 
if( i .eg. idiv ) thdx = ( 1.-sdispl)*pi 
else 
thax = todx 
endif 
THDX = DTHXI(I) 
thdxd = (180/pi)*thax 
write(6,*; i, thdxd 
DAW(I) = abs ( SIN( THDX ) } 
DAWS = DXKS + DXW(1) 
DO 20 I = 1,IDIV 
DR(I) = DXB(I) * AOL / DAWS 
v 


Jo 


_ 
OOO OW O&O oe 08'08' Ge 00 0 .On'00: 
NOV WN — ODO GNU Binks 


} 
i 


oe ee ei ane! ae ee oe oe: 


DxXI = 0. 

DO 3 I= 1, IDI 

DRI = DXI + DxX(1) 

XX = DXI + AOST 

CALL TAINT(A0,S0,XX,YY,N,3,NER,MON) 
XP(I) = XX * XX - YY * YY + XSING 
YP(I) = 2. * XX * YY + YSING 

RETURN 

END 









SOURCE PROGRAM 
airfgr. 
= 





ie woe ee 
















SOURCE TEXT 


SUBROUTINE TAINT(XTAB, PTAB,X%,FX,N,E,NER,MON) 


Crserwasnsseesansteeeesnetereesanecerres RACeeRaateseeaanaserseeaatrese 
rare * 
ce SUBROUTINE TAINT * 
ee . 
CEAaeeeaeesae eer een es ee eee sRAeSHEIDALSSHSSRERUSESSRADE DEES A ERA SET EAE 
PARAMETER (1X=300, KX2100) 

DIMENSION XTAB(1),FTAB(1),7(10) ,C(10) 








MASA - AMES SUBROUTINE FOR POLYNOMIAL INTERPOLATION 
OF A TABULATED FUNCTION 


aaaa 


IPF(N-K) 1,1 , 2 
NER = 2 
RETURN 


pos rs 1, NM 
T¥(XTAB(I) - XTABC(I+1)}) 9,22,10 


11 NER = 3 
RETURN 
93= 3-1 
GO TO 8 
i 10 J = J+1 
7239} 8 CONTINUE 
“2407 MON = 1 
T241 IF(J) 12,6, 6 
242] 12 MON = 2 
243 7poiste1,N 
244 IF(X - 2LTAB(I)) 14,14,13 
7245; wget 
7246} Go To 18 
2475 13 CONTINUE 
_248) GO To 15 
249} 6polérel,N 
250 | IF(X-2TAB(1)) 16,17,17 
p25} wwIJ=1 
282! GO To 18 : 
253! 16 CONTINUE 
7254; iSJen 
255° 18 J = J - (K-1) / 2 
. 256: IF(J) 19,19,20 
257 i 9Je1 
258 : 2°Me J+ K 
“289 | IE(M - N) 21,21,22 
260: 223=3-1 
261: GO TO 20 
262 | 21 KP1 = kX +1 
"263 | JSAVE = J 
264 | 26 pO 23 L = 1, KP 
__ 265 : C(L) = X - XTAB(J) 
_266 | T(L) = FTAB()) 
_267 | 23.3 = J+1 
268 ; DO 24.3 = 1,K 
269 | 1 J+1 
270 | 25 TT) = (CEZ) TET -C( I) 9 TI) 7 (C(T)-C(L)) 
271; I= Il 
_ 272 | IF(I-FP1) 25,25,24 
_ 273 | 24 CONTINCE 
274; FX = T(KP1) 
275, NER = 1 
276 RETURN 


277 | END 











SOURCE PROGRAM DATE —«- 1/20/90 
airfgr. | TME 11:45:57 am | 


SUBROUTINE SING(N2,N,X,2,X%LE ,YLE, TEA, TES, XSING , YSING,CBD) 
CeCOOCAR RES S CER NHAC EH ERENSSERERRSTTOTRARESCERE TESST SRE ERAS ET ENARIES 
ct * 
c SUBROUTINE SING : ij 
ce * 
CleReneeaeer aN eed eERASOOSERAKHT TE RTERASOOURRECCOSERARSSOTRENTHO ER ERS 


PARAMETER (1X*300,KX#100) 


THIS SUBROUTINE COMPUTES SINGULAR POINT LOCATIONS. 


DIMENSION X(2) , 2(2) 

NU = N2 

Nl = N24 2 

N3 N2-1 

zu (N21) 

zu Z(N1) 

12 X(N2) 

22 282) 

x3 2(N3) 

23 2(N3) 

D1 420°* 2- 21 

D2 m2 «8 2- 22 

D3 42 - x1 

mo 22 - 21 

D5 x3 ** 2 - x2 

ps 23 ** 2 - 21 

D7 13 -- x1 

os 23 - 21 

B= (D7 * ( Di + D2) - D3*(DS+D6))/(2.*(D7*D4-D3*DB) ) 
IP(ABS(D3).LT.ABS(D7)) GO TO 10 
A= (D1 + D2 - 2. * B* D4) / (2. * D3) 
GO TO 20 

A= (D5 + D6 - 2. * B * D8) /(¢ 2. * D7) 
CONTINUE 

R = SORT((X2-A)** 2 + (Z2-B)**2) 
RLE = X(NU) 

YLE = 2(NU) 

CHD = X(1) - X(NU) 

AZ = (2(2)-2(1)) /(K(2) - XC1)) 
Al = (2(N)-2(N-1))/(%(N)-X(N-1)) 
TES = .5 * (Al + A2) 

TEA = A2 - Al 

TEA = TEA * 57.29578 

USING = (A+XLE) /2. 

YSING * (B+YLE) / 2. 

RETURN 

END 


L 


eeeeae 


oo 
n= 


al ada 


ie A PEA SEES EN BEES Sy SE a 


VBWN=—CUONAV AWN —O. 
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SOURCE PROCRAM 
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SUBROUTINE AIRFOL( AORAT , adexp ,sdispl) 


| CH OOAeENADS OEE ATRESSEERASHORAEERESSHRETKHSSCOR ERAS SOLER ERHSCEUREEESOOE 


ce SUBROUTINE ATRFOL e 
Ces eeccaatateeansrat sean ates eesaeetseesestettesanttitenerisetraastasese 
i PARAMETER (1X=300,KX=100) 

COMMON/GRID1/X( IE, KX),2(1IZ,KX) 

COMMON/DGRID/DT , IMAX, KMAZ, ITEL, ITEU 

COMMON /YSYM/ISYM 

DIMENSION SO(IX,4),AQ(IX,4),YO(IX,4),XP(IX),YP(13), 

LE(IX) FCI), BOC IX) 

DIMENSION XL(IZ), XUCIZ), YL(IR), YUCIX), 

2 AX(IX), Y¥(1X) 


DATA (B0( 2), T=1,32)/2. ,2.0424,1. 0836 ,1.2270,1.1715,2.2175,1.2651, 
21.3145,1.3659,1.4199,2 .4755,1.5349,1.5973,1.6636,1.7342,1.8099, 
23, 85246,2.9799,2.0764,2.2823,2.3022,2.4342,2. 5653.2. 7597,2. 9546, 
33.2106,3.5141,3.9019,4.4219,5.1687,8.3632,8.6809/ 


READ(S,*) ISYM, IBMAX 
IF( ISYM .£Q. 1 ) THEN 
DO 101 1*1,IBMAX 
READ(S,*) XU(1), YU(1) 
ELSE 

Do 102 I=1,1BMAX 
READ(S,*) XU(1), YU(I), YL(I) 
CONTINUE 

ENDIF 

IF( ISYM .£Q. 1 ) TBEN 
DO 103 1*#1,IBMAX 
YL(I) = = YUCT) 

ENDIF 


DO 1000 I=), IBMAX 
IU © I + IBMAX 
AX(IV) © XLT) 
YY(IU) = YU(I) 

IL © IBMAX - I+ 1 
X(T) XU TL) 
Y¥(1) = YL(IL) 
CONTINUE 


IBMAX2 = 2*IBMAX 
DO 1010 = T=1, IBMAX2 
AP(I) = XX(1) 

YP(I) = YY(1) 
CONTINUE 

FNU = IBMAX 

FNL = IBMAX 


TRIS SUBROVTINE GENERATES SKEAR PARABOLIC C-GRID 

THE FOLLOWING SUBROUTINES ARE RELATED TO THE GRID GENERATION 
WRAP SING 
TABINT CLUSTR 
TAINT STRTC 


aanaaana 


AQ(1,1) = 1. 

pO ®@ I= 2, IBMAX 
AO(T,2) = ( AQ(I~1,1) * AORAT )**AQEXP 
PI = 4. * ATAN(1.) 

NU = FRU 

NL © FNL 

No* NU + NL 

IT = ITEL 

IE = YTEU 

ILE= ( ITEL + ITEU ) / 2 
II = 

Ja 

11P1 

TIM1 

11JJ 

IIJJ2 © 11 * ( JJ-2) 


SCALE = 1./ ( XP(1) - XP(NL) ) 

DO 3333 1*1,N 

XP(1) = XP(I) * SCALE 

YP(I) = YP(I) * SCALE 

CONTINUE 

CALL SING(NU,N,XP,YP,XLE,ZLE, TEA, TES, XSING, YSING, CED) 
CALL TABINT(XP,YP,XSING, YSING,N,8d15p1) 
NBODY © IE + 1 - IT 

DO 6791 I = 1 , NBODY 

Ler-1 

E(IT+L) = XP(1) 

F(IT+L) = YP(I) 

IEP1 = IE + 1 

SLOPT «© TES * .75 

po 436 1 = JEP] , II 

Ile I 41 - JE 

E(1) = AQ(11,1) 

DXI= 1. / 48. 

D = 4. / 3. * (E(J) - .25) 

F(I) © FCIE) + SLOPT * ALOG(D) / D 

L= 11Pl - 1 

E(L) © E(1) 

F(L) = F(IT) + SLOPT © ALOG(D)/D 

WRITE (6,439) 

FORMAT(22,3B 1,19X,1NK,29X, 187) 

WRITE (6,37) (1,B(1),8(2),1 = 1, 27) 
CALL WRAP(II,J3J,XSING,YSING,E,F,$0,A0,Y0) 
po 10 J* 2, 33 

poior*i1, 11 

ACL,J~2) © AO(T,1)¢°2 - (S0(1,1)+¥0(9,1))**2 
Z(T,I-1) © 2. © AO(T,1) * ($0(1,1)4¥0(3,1)) 
RETURN 

FORMAT(I5,2F20.8) 

END 





SOURCE PROGRAM 
airfgr.f 


ee: 


SUBROUTINE CLUSTR(DS) 
CHSC HE RAAHSSEARASHHSARAAHH SHUNT TLESSARAEKECSSERAEHOSSCENEREDECEEAAAOOE 
.@. 
SUBROUTINE CLUSTR ‘ 
* 
37 Cee RAASS SHAR KHSSSHRAAFSSSHAATAHSSSSAEAAASSSHRE RASCH ENEAASSOSHUREROCCEREE 
8) PARAMETER (1X=300,KX=100) 
COMMON/GRID1/X (1X, KX}, 2(1%, KX) 
COMMON/DGRID/DT, IMAX , KMAX, I TEL, ITEU 
DIMENSION S(IX),XP(IX),YP(IX),R(1X) 


FAIS SUBROUTINE CLUSTERS A GIVER X,%2 GRID SUCHE THAT THE FIRST POINT IS AT 


DO 100 1 © 1, IMAX 
sq) = 0. 
EP(1) = X(1,1) 
YP(1l) © 2(2,2) 
DO 10 R= 2, KMAX 
ZP(K) = X(1,K) 
YP(K) © 2(1,K) 
10 $(K) = SORT( (XP(K)-XP(R-1))**2+(YP(R)-YP(R-1))**2) 
148(K-1) 
SUMDX = S(KMAX) 
CALL STRTCH(SUMDX ,DS,F1,KMAX , FACTOR) 
WRITE (6,200) 1, FACTOR 
R(1) = 0. 
DR = DS 
po 20 K = 2 , KMAX 
R(R) = R(K-1) + DR 
DR = DR * FACTOR 
CONTINUE 
RLAST « 1. / R(KMAX) 
DO 30 K = 2 , KMAX 
Rl = R(K) * RLAST © S(KMAX) 
CALL TAINT(S,XP,R1,XP1,RMAX,3,NER,MON) 
X(I,K) = XPL 
CALL TAINT(S,YP,R1,Y22, KMAX,3,NER,MON) 
Z(1,R) * YP1 
CONTINUE 
CONTINUE 
WRITE (6,115) 
DO 110 1 = 2, IMAX 
DX = X(1,2) - X(T,1) 
DY = 2(1,2) - 2(1,1) 
DN = SORT/DX*DX+DY *DY) 
WRITE(6,120) I , DX , DY , DN 
110 CONTINUE 
RETURN 
115 FORMAT(SX,6HNORMAL,1X,8HDISTANCE,3B AT,4B THE,5H WALL,/ 
1,5B 1,8X,2BDX, 8X, 2HDY , 8X, 2HDN,//) 
120 FORMAT(I5,3F10.5) 
200 FORMAT(I5,F10.5) 
END 










485 

87 

ie a 
489"! 











ce 
ce 
ct 


SOURCE PROCRAM 


: 
airfgr.f [TIME 11:45:57 am | 


TIME §=14:45:57 am 












SUBROUTINE STRTCB({SUMDX , DX), F1,N1,R) 


4a Cooeeeennnccerens ceerewenSeeeeenn see SeERKSSCCERN AAS SORE ECAESESPRAAEOT 


SUBROUTINE STRITCH re 
* 


CHUNaeoeekaansese cane eseSA RENE OERERRASS SORENESS CERERESOORTESSS SO OEEE 


10 


PARAMETER (1%=300,KX=100) 


THIS SUBROUTINE DETERMINES A GEOMETRIC . 
PROGRESSION OF GRID SPACING BETWEEN } AND Nl SUH TSAT 
BSUMSDX) SQUALS SUMDX. THE HATIO BETWEEN SUCCESSIVE 
SPACINGS IS R. 

NeNl- 1 

R-1.5 

El = 1.E~04 

E2 = 1.E-04 

po 10 L=1, 50 

Fe (R-1) * SUMDX - DX1*(R**N-1) 

FP = SUMDI - DXl * FLOAT(N) * R ** (N-1) 

RITER = R - F/ FP 


“2F(L.E-02. 17. RITER.AND.RITER.LT.1.} RITER = 1. 


FHL. .LT.RITER.AND.RITER-LT.100.) RITER=.01 
IF(ABS(R-RITER).LT. R*E1) GO To 1 

R = RITER 

CONTINUE 

R = 1.0001 

OX2 = DUTOT/FLOAT(N1-1} 

RETURN 

R= RITER 

RETURN 

END 





SOURCE PROCRAM 


airfgr.f 


ROTGRID({ IMAZ, KMAZ, DALFA) 
CASHRAERATSSHSRATASSSERAARSESORARATHSSERALHESOERTLASSEERRASHCEERLERGES 
i e 
SUBROUTINE ROTGID ne * 
t 
Clee KSSeeaA RH SeeleRkaASSCORENAASSCEREKSSSORAAAHLOTHLAKTHSSHAETERSOERERE 
2) PARAMETER (IX=300,K2=100) 
ic ROTATE GRID IN THE CLOCKWISE DIRECTION BY AN AMOUNT DALFA 
COMMON/GRID1/X( IX, KX), 2(IX,KX) 
CA = COS(DALFA) 
SA = - SIN(DALFA) 
DO 20 K = 1, KMAX 
DO 201 = 2, IMAX 
Xl = x(1,2) 
Zl = 2(I,K) 
(I,K) = Xl * CA - Zi * SA 
Z(1,K) = 21 ¢ CA + x1 * SA 
RETURN 
END 
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